matrix.hpp

Go to the documentation of this file.

//============================================================================
// MCKL/include/mckl/core/matrix.hpp
//----------------------------------------------------------------------------
// MCKL: Monte Carlo Kernel Library
//----------------------------------------------------------------------------
// Copyright (c) 2013-2018, Yan Zhou
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
//   Redistributions of source code must retain the above copyright notice,
//   this list of conditions and the following disclaimer.
//
//   Redistributions in binary form must reproduce the above copyright notice,
//   this list of conditions and the following disclaimer in the documentation
//   and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//============================================================================

#ifndef MCKL_CORE_MATRIX_HPP
#define MCKL_CORE_MATRIX_HPP

#include <mckl/internal/common.hpp>
#include <mckl/core/is_equal.hpp>
#include <mckl/core/iterator.hpp>

namespace mckl {

template <typename T, MatrixLayout Layout = ColMajor,
    typename Alloc = Allocator<T>>
class Matrix
{
    using layout_dispatch = std::integral_constant<MatrixLayout, Layout>;
    using row_major = std::integral_constant<MatrixLayout, RowMajor>;
    using col_major = std::integral_constant<MatrixLayout, ColMajor>;

    Vector<T, Alloc> data_;

  public:
    using value_type = T;
    using allocator_type = Alloc;
    using size_type = std::size_t;
    using difference_type = std::ptrdiff_t;
    using reference = value_type &;
    using const_reference = const value_type &;
    using pointer = value_type *;
    using const_pointer = const value_type *;

    using iterator = pointer;
    using const_iterator = const_pointer;
    using reverse_iterator = std::reverse_iterator<iterator>;
    using const_reverse_iterator = std::reverse_iterator<const_iterator>;

    using row_iterator =
        std::conditional_t<Layout == RowMajor, pointer, StepIterator<pointer>>;
    using const_row_iterator = std::conditional_t<Layout == RowMajor,
        const_pointer, StepIterator<const_pointer>>;
    using reverse_row_iterator = std::reverse_iterator<row_iterator>;
    using const_reverse_row_iterator =
        std::reverse_iterator<const_row_iterator>;

    using col_iterator =
        std::conditional_t<Layout == ColMajor, pointer, StepIterator<pointer>>;
    using const_col_iterator = std::conditional_t<Layout == ColMajor,
        const_pointer, StepIterator<const_pointer>>;
    using reverse_col_iterator = std::reverse_iterator<col_iterator>;
    using const_reverse_col_iterator =
        std::reverse_iterator<const_col_iterator>;

    using row_range = Range<row_iterator>;
    using const_row_range = Range<const_row_iterator>;
    using reverse_row_range = Range<reverse_row_iterator>;
    using const_reverse_row_range = Range<const_reverse_row_iterator>;

    using col_range = Range<col_iterator>;
    using const_col_range = Range<const_col_iterator>;
    using reverse_col_range = Range<reverse_col_iterator>;
    using const_reverse_col_range = Range<const_reverse_col_iterator>;

    using transpose_type =
        Matrix<T, Layout == RowMajor ? ColMajor : RowMajor, Alloc>;

    Matrix() noexcept(
        std::is_nothrow_default_constructible<Vector<T, Alloc>>::value)
        : nrow_(0), ncol_(0)
    {
    }

    Matrix(size_type nrow, size_type ncol)
        : data_(nrow * ncol), nrow_(nrow), ncol_(ncol)
    {
    }

    Matrix(const Matrix &) = default;

    Matrix(Matrix &&other) noexcept(
        std::is_nothrow_move_constructible<Vector<T, Alloc>>::value)
        : data_(std::move(other.data_)), nrow_(other.nrow_), ncol_(other.ncol_)
    {
        other.nrow_ = 0;
        other.ncol_ = 0;
    }

    Matrix &operator=(const Matrix &) = default;

    Matrix &operator=(Matrix &&other) noexcept(
        noexcept(data_.swap(other.data_)))
    {
        if (this != &other) {
            data_.swap(other.data_);
            std::swap(nrow_, other.nrow_);
            std::swap(ncol_, other.ncol_);
        }

        return *this;
    }

    explicit operator transpose_type() const
    {
        transpose_type mat(nrow_, ncol_);

        if (Layout == RowMajor) {
            for (size_type i = 0; i != nrow_; ++i) {
                for (size_type j = 0; j != ncol_; ++j) {
                    mat(i, j) = operator()(i, j);
                }
            }
        }
        if (Layout == ColMajor) {
            for (size_type j = 0; j != ncol_; ++j) {
                for (size_type i = 0; i != nrow_; ++i) {
                    mat(i, j) = operator()(i, j);
                }
            }
        }

        return mat;
    }

    allocator_type get_allocator() const { return data_.get_allocator(); }

    iterator begin() { return data(); }

    const_iterator begin() const { return cbegin(); }

    const_iterator cbegin() const { return data(); }

    iterator end() { return begin() + nrow_ * ncol_; }

    const_iterator end() const { return cend(); }

    const_iterator cend() const { return begin() + nrow_ * ncol_; }

    reverse_iterator rbegin() { return reverse_iterator(end()); }

    const_reverse_iterator rbegin() const { return crbegin(); }

    const_reverse_iterator crbegin() const
    {
        return const_reverse_iterator(cend());
    }

    reverse_iterator rend() { return reverse_iterator(begin()); }

    const_reverse_iterator rend() const { return crend(); }

    const_reverse_iterator crend() const
    {
        return const_reverse_iterator(cbegin());
    }

    row_iterator row_begin(size_type i)
    {
        return row_begin_dispatch(i, layout_dispatch());
    }

    const_row_iterator row_begin(size_type i) const { return row_cbegin(i); }

    const_row_iterator row_cbegin(size_type i) const
    {
        return row_begin_dispatch(i, layout_dispatch());
    }

    row_iterator row_end(size_type i) { return row_begin(i) + ncol_; }

    const_row_iterator row_end(size_type i) const { return row_cend(i); }

    const_row_iterator row_cend(size_type i) const
    {
        return row_cbegin(i) + ncol_;
    }

    reverse_row_iterator row_rbegin(size_type i)
    {
        return reverse_row_iterator(row_end(i));
    }

    const_reverse_row_iterator row_rbegin(size_type i) const
    {
        return row_crbegin(i);
    }

    const_reverse_row_iterator row_crbegin(size_type i) const
    {
        return const_reverse_row_iterator(row_cend(i));
    }

    reverse_row_iterator row_rend(size_type i)
    {
        return reverse_row_iterator(row_begin(i));
    }

    const_reverse_row_iterator row_rend(size_type i) const
    {
        return row_crend(i);
    }

    const_reverse_row_iterator row_crend(size_type i) const
    {
        return const_reverse_row_iterator(row_cbegin(i));
    }

    col_iterator col_begin(size_type i)
    {
        return col_begin_dispatch(i, layout_dispatch());
    }

    const_col_iterator col_begin(size_type j) const { return col_cbegin(j); }

    const_col_iterator col_cbegin(size_type j) const
    {
        return col_begin_dispatch(j, layout_dispatch());
    }

    col_iterator col_end(size_type j) { return col_begin(j) + nrow_; }

    const_col_iterator col_end(size_type j) const { return col_cend(j); }

    const_col_iterator col_cend(size_type j) const
    {
        return col_cbegin(j) + nrow_;
    }

    reverse_col_iterator col_rbegin(size_type j)
    {
        return reverse_col_iterator(col_end(j));
    }

    const_reverse_col_iterator col_rbegin(size_type j) const
    {
        return col_crbegin(j);
    }

    const_reverse_col_iterator col_crbegin(size_type j) const
    {
        return const_reverse_col_iterator(col_cend(j));
    }

    reverse_col_iterator col_rend(size_type j)
    {
        return reverse_col_iterator(col_begin(j));
    }

    const_reverse_col_iterator col_rend(size_type j) const
    {
        return col_crend(j);
    }

    const_reverse_col_iterator col_crend(size_type j) const
    {
        return const_reverse_col_iterator(col_cbegin(j));
    }

    row_range row(size_type i) { return row_range(row_begin(i), row_end(i)); }

    const_row_range row(size_type i) const
    {
        return const_row_range(row_begin(i), row_end(i));
    }

    const_row_range crow(size_type i) const
    {
        return const_row_range(row_cbegin(i), row_cend(i));
    }

    reverse_row_range rrow(size_type i)
    {
        return reverse_row_range(row_rbegin(i), row_rend(i));
    }

    const_reverse_row_range rrow(size_type i) const
    {
        return const_reverse_row_range(row_rbegin(i), row_rend(i));
    }

    const_reverse_row_range crrow(size_type i) const
    {
        return const_reverse_row_range(row_crbegin(i), row_crend(i));
    }

    col_range col(size_type i) { return col_range(col_begin(i), col_end(i)); }

    const_col_range col(size_type i) const
    {
        return const_col_range(col_begin(i), col_end(i));
    }

    const_col_range ccol(size_type i) const
    {
        return const_col_range(col_cbegin(i), col_cend(i));
    }

    reverse_col_range rcol(size_type i)
    {
        return reverse_col_range(col_rbegin(i), col_rend(i));
    }

    const_reverse_col_range rcol(size_type i) const
    {
        return const_reverse_col_range(col_rbegin(i), col_rend(i));
    }

    const_reverse_col_range crcol(size_type i) const
    {
        return const_reverse_col_range(col_crbegin(i), col_crend(i));
    }

    reference at(size_type i, size_type j)
    {
        runtime_assert<std::out_of_range>(
            i < nrow_, "**Matrix::at** row index out of range");
        runtime_assert<std::out_of_range>(
            j < ncol_, "**Matrix::at** column index out of range");

        return at_dispatch(i, j, layout_dispatch());
    }

    const_reference at(size_type i, size_type j) const
    {
        runtime_assert<std::out_of_range>(
            i < nrow_, "**Matrix::at** row index out of range");
        runtime_assert<std::out_of_range>(
            j < ncol_, "**Matrix::at** column index out of range");

        return at_dispatch(i, j, layout_dispatch());
    }

    reference operator()(size_type i, size_type j)
    {
        return at_dispatch(i, j, layout_dispatch());
    }

    const_reference operator()(size_type i, size_type j) const
    {
        return at_dispatch(i, j, layout_dispatch());
    }

    pointer data() { return data_.data(); }

    const_pointer data() const { return data_.data(); }

    pointer row_data(size_type i)
    {
        return row_data_dispatch(i, layout_dispatch());
    }

    const_pointer row_data(size_type i) const
    {
        return row_data_dispatch(i, layout_dispatch());
    }

    size_type row_stride() const
    {
        return row_stride_dispatch(layout_dispatch());
    }

    pointer col_data(size_type j)
    {
        return col_data_dispatch(j, layout_dispatch());
    }

    const_pointer col_data(size_type j) const
    {
        return col_data_dispatch(j, layout_dispatch());
    }

    size_type col_stride() const
    {
        return col_stride_dispatch(layout_dispatch());
    }

    static constexpr MatrixLayout layout() { return Layout; }

    bool empty() const { return nrow_ == 0 || ncol_ == 0; }

    size_type nrow() const { return nrow_; }

    size_type ncol() const { return ncol_; }

    size_type ldim() const { return ldim_dispatch(layout_dispatch()); }

    size_type size() const { return nrow_ * ncol_; }

    void reserve(size_type n, size_type m) { data_.reserve(n * m); }

    void resize(size_type nrow, size_type ncol)
    {
        if (nrow * ncol == 0) {
            data_.clear();
            nrow_ = nrow;
            ncol_ = ncol;
            return;
        }

        if (nrow == nrow_ && ncol_ == ncol) {
            return;
        }

        resize_dispatch(nrow, ncol, layout_dispatch());
    }

    void shrink_to_fit() { data_.shrink_to_fit(); }

    void clear()
    {
        data_.clear();
        nrow_ = 0;
        ncol_ = 0;
    }

    template <typename InputIter>
    void push_back_row(InputIter first)
    {
        resize(nrow_ + 1, ncol_);
        std::copy_n(first, ncol_, row_begin(nrow_ - 1));
    }

    template <typename InputIter>
    void push_back_col(InputIter first)
    {
        resize(nrow_, ncol_ + 1);
        std::copy_n(first, nrow_, col_begin(ncol_ - 1));
    }

    void swap(Matrix &other) noexcept(noexcept(data_.swap(other.data_)))
    {
        std::swap(nrow_, other.nrow_);
        std::swap(ncol_, other.ncol_);
        data_.swap(other.data_);
    }

    friend void swap(Matrix &m1, Matrix &m2) noexcept(noexcept(m1.swap(m2)))
    {
        m1.swap(m2);
    }

    friend bool operator==(const Matrix &m1, const Matrix &m2)
    {
        if (m1.nrow_ != m2.nrow_) {
            return false;
        }
        if (m1.ncol_ != m2.ncol_) {
            return false;
        }
        return m1.data_ == m2.data_;
    }

    friend bool operator!=(const Matrix &m1, const Matrix &m2)
    {
        return !(m1 == m2);
    }

    friend bool is_equal(const Matrix &m1, const Matrix &m2)
    {
        if (m1.nrow_ != m2.nrow_) {
            return false;
        }
        if (m1.ncol_ != m2.ncol_) {
            return false;
        }
        return is_equal(m1.data_, m2.data_);
    }

  private:
    std::size_t nrow_;
    std::size_t ncol_;

    // Layout == RowMajor

    size_type ldim_dispatch(row_major) const { return ncol_; }

    pointer row_begin_dispatch(size_type i, row_major) { return row_data(i); }

    const_pointer row_begin_dispatch(size_type i, row_major) const
    {
        return row_data(i);
    }

    StepIterator<pointer> col_begin_dispatch(size_type j, row_major)
    {
        return StepIterator<pointer>(
            col_data(j), static_cast<difference_type>(col_stride()));
    }

    StepIterator<const_pointer> col_begin_dispatch(
        size_type j, row_major) const
    {
        return StepIterator<const_pointer>(
            col_data(j), static_cast<difference_type>(col_stride()));
    }

    reference at_dispatch(size_type i, size_type j, row_major)
    {
        return data_[i * ncol_ + j];
    }

    const_reference at_dispatch(size_type i, size_type j, row_major) const
    {
        return data_[i * ncol_ + j];
    }

    size_type row_stride_dispatch(row_major) const { return 1; }

    pointer row_data_dispatch(size_type i, row_major)
    {
        return data() + i * ncol_;
    }

    const_pointer row_data_dispatch(size_type i, row_major) const
    {
        return data() + i * ncol_;
    }

    size_type col_stride_dispatch(row_major) const { return ncol_; }

    pointer col_data_dispatch(size_type j, row_major) { return data() + j; }

    const_pointer col_data_dispatch(size_type j, row_major) const
    {
        return data() + j;
    }

    void resize_dispatch(size_type nrow, size_type ncol, row_major)
    {
        if (ncol == ncol_) {
            data_.resize(nrow * ncol);
            nrow_ = nrow;
            return;
        }

        Matrix tmp(nrow, ncol);
        const size_type n = std::min(nrow, nrow_);
        const size_type m = std::min(ncol, ncol_);
        for (size_type i = 0; i != n; ++i) {
            std::copy_n(row_data(i), m, tmp.row_data(i));
        }
        swap(tmp);
    }

    // Layout == ColMajor

    size_type ldim_dispatch(col_major) const { return nrow_; }

    StepIterator<pointer> row_begin_dispatch(size_type i, col_major)
    {
        return StepIterator<pointer>(
            row_data(i), static_cast<difference_type>(row_stride()));
    }

    StepIterator<const_pointer> row_begin_dispatch(
        size_type i, col_major) const
    {
        return StepIterator<const_pointer>(
            row_data(i), static_cast<difference_type>(row_stride()));
    }

    pointer col_begin_dispatch(size_type j, col_major) { return col_data(j); }

    const_pointer col_begin_dispatch(size_type j, col_major) const
    {
        return col_data(j);
    }

    reference at_dispatch(size_type i, size_type j, col_major)
    {
        return data_[j * nrow_ + i];
    }

    const_reference at_dispatch(size_type i, size_type j, col_major) const
    {
        return data_[j * nrow_ + i];
    }

    size_type row_stride_dispatch(col_major) const { return nrow_; }

    pointer row_data_dispatch(size_type i, col_major) { return data() + i; }

    const_pointer row_data_dispatch(size_type i, col_major) const
    {
        return data() + i;
    }

    size_type col_stride_dispatch(col_major) const { return 1; }

    pointer col_data_dispatch(size_type j, col_major)
    {
        return data() + j * nrow_;
    }

    const_pointer col_data_dispatch(size_type j, col_major) const
    {
        return data() + j * nrow_;
    }

    void resize_dispatch(size_type nrow, size_type ncol, col_major)
    {
        if (nrow == nrow_) {
            data_.resize(nrow * ncol);
            ncol_ = ncol;
            return;
        }

        Matrix tmp(nrow, ncol);
        const size_type n = std::min(nrow, nrow_);
        const size_type m = std::min(ncol, ncol_);
        for (size_type i = 0; i != m; ++i) {
            std::copy_n(col_data(i), n, tmp.col_data(i));
        }
        swap(tmp);
    }
}; // class Matrix

template <typename CharT, typename Traits, typename T, MatrixLayout Layout,
    typename Alloc>
inline std::basic_ostream<CharT, Traits> &operator<<(
    std::basic_ostream<CharT, Traits> &os, const Matrix<T, Layout, Alloc> &mat)
{
    if (!os) {
        return os;
    }

    os << mat.nrow() << ' ' << mat.ncol();
    if (!os) {
        return os;
    }

    for (std::size_t j = 0; j != mat.ncol(); ++j) {
        for (std::size_t i = 0; i != mat.nrow(); ++i) {
            os << ' ' << mat(i, j);
        }
    }

    return os;
}

template <typename CharT, typename Traits, typename T, MatrixLayout Layout,
    typename Alloc>
inline std::basic_istream<CharT, Traits> &operator>>(
    std::basic_istream<CharT, Traits> &is, Matrix<T, Layout, Alloc> &mat)
{
    if (!is) {
        return is;
    }

    std::size_t nrow = 0;
    std::size_t ncol = 0;
    is >> nrow >> std::ws >> ncol;
    if (!is) {
        return is;
    }

    Matrix<T, Layout, Alloc> tmp(nrow, ncol);
    for (std::size_t j = 0; j != mat.ncol(); ++j) {
        for (std::size_t i = 0; i != mat.nrow(); ++i) {
            is >> std::ws >> tmp(i, j);
        }
    }
    if (is) {
        mat = std::move(tmp);
    }

    return is;
}

} // namespace mckl

#endif // MCKL_CORE_MATRIX_HPP