matrix.h

Go to the documentation of this file.

/***************************************************************************
 *  include/stxxl/bits/containers/matrix.h
 *
 *  Part of the STXXL. See http://stxxl.sourceforge.net
 *
 *  Copyright (C) 2010-2011 Raoul Steffen <[email protected]>
 *
 *  Distributed under the Boost Software License, Version 1.0.
 *  (See accompanying file LICENSE_1_0.txt or copy at
 *  http://www.boost.org/LICENSE_1_0.txt)
 **************************************************************************/

#ifndef STXXL_CONTAINERS_MATRIX_HEADER
#define STXXL_CONTAINERS_MATRIX_HEADER

#include <stxxl/bits/containers/vector.h>
#include <stxxl/bits/common/counting_ptr.h>
#include <stxxl/bits/mng/block_scheduler.h>
#include <stxxl/bits/containers/matrix_arithmetic.h>

STXXL_BEGIN_NAMESPACE

//! \defgroup matrix matrix
//! Efficient external memory matrix operations
//! \ingroup stlcont
//! \{

/* index-variable naming convention:
 * [MODIFIER_][UNIT_]DIMENSION[_in_[MODIFIER_]ENVIRONMENT]
 *
 * e.g.:
 * block_row = number of row measured in rows consisting of blocks
 * element_row_in_block = number of row measured in rows consisting of elements in the (row of) block(s)
 *
 * size-variable naming convention:
 * [MODIFIER_][ENVIRONMENT_]DIMENSION[_in_UNITs]
 *
 * e.g.
 * height_in_blocks
 */

// forward declaration
template <typename ValueType, unsigned BlockSideLength>
class matrix;

//! external column-vector container for matrix multiplication
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
template <typename ValueType>
class column_vector : public vector<ValueType>
{
public:
    typedef vector<ValueType> vector_type;
    typedef typename vector_type::size_type size_type;

    using vector_type::size;

    //! \param n number of elements
    column_vector(size_type n = 0)
        : vector_type(n) { }

    column_vector operator + (const column_vector& right) const
    {
        assert(size() == right.size());
        column_vector res(size());
        for (size_type i = 0; i < size(); ++i)
            res[i] = (*this)[i] + right[i];
        return res;
    }

    column_vector operator - (const column_vector& right) const
    {
        assert(size() == right.size());
        column_vector res(size());
        for (size_type i = 0; i < size(); ++i)
            res[i] = (*this)[i] - right[i];
        return res;
    }

    column_vector operator * (const ValueType scalar) const
    {
        column_vector res(size());
        for (size_type i = 0; i < size(); ++i)
            res[i] = (*this)[i] * scalar;
        return res;
    }

    column_vector& operator += (const column_vector& right)
    {
        assert(size() == right.size());
        for (size_type i = 0; i < size(); ++i)
            (*this)[i] += right[i];
        return *this;
    }

    column_vector& operator -= (const column_vector& right)
    {
        assert(size() == right.size());
        for (size_type i = 0; i < size(); ++i)
            (*this)[i] -= right[i];
        return *this;
    }

    column_vector& operator *= (const ValueType scalar)
    {
        for (size_type i = 0; i < size(); ++i)
            (*this)[i] *= scalar;
        return *this;
    }

    void set_zero()
    {
        for (typename vector_type::iterator it = vector_type::begin(); it != vector_type::end(); ++it)
            *it = 0;
    }
};

//! external row-vector container for matrix multiplication
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
template <typename ValueType>
class row_vector : public vector<ValueType>
{
public:
    typedef vector<ValueType> vector_type;
    typedef typename vector_type::size_type size_type;

    using vector_type::size;

    //! \param n number of elements
    row_vector(size_type n = 0)
        : vector_type(n) { }

    row_vector operator + (const row_vector& right) const
    {
        assert(size() == right.size());
        row_vector res(size());
        for (size_type i = 0; i < size(); ++i)
            res[i] = (*this)[i] + right[i];
        return res;
    }

    row_vector operator - (const row_vector& right) const
    {
        assert(size() == right.size());
        row_vector res(size());
        for (size_type i = 0; i < size(); ++i)
            res[i] = (*this)[i] - right[i];
        return res;
    }

    row_vector operator * (const ValueType scalar) const
    {
        row_vector res(size());
        for (size_type i = 0; i < size(); ++i)
            res[i] = (*this)[i] * scalar;
        return res;
    }

    template <unsigned BlockSideLength>
    row_vector operator * (const matrix<ValueType, BlockSideLength>& right) const
    { return right.multiply_from_left(*this); }

    ValueType operator * (const column_vector<ValueType>& right) const
    {
        ValueType res = 0;
        for (size_type i = 0; i < size(); ++i)
            res += (*this)[i] * right[i];
        return res;
    }

    row_vector& operator += (const row_vector& right)
    {
        assert(size() == right.size());
        for (size_type i = 0; i < size(); ++i)
            (*this)[i] += right[i];
        return *this;
    }

    row_vector& operator -= (const row_vector& right)
    {
        assert(size() == right.size());
        for (size_type i = 0; i < size(); ++i)
            (*this)[i] -= right[i];
        return *this;
    }

    row_vector& operator *= (const ValueType scalar)
    {
        for (size_type i = 0; i < size(); ++i)
            (*this)[i] *= scalar;
        return *this;
    }

    void set_zero()
    {
        for (typename vector_type::iterator it = vector_type::begin(); it != vector_type::end(); ++it)
            *it = 0;
    }
};

//! Specialized swappable_block that interprets uninitialized as containing zeros.
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
//! \tparam BlockSideLength side length of a matrix block
//!
//! When initializing, all values are set to zero.
template <typename ValueType, unsigned BlockSideLength>
class matrix_swappable_block : public swappable_block<ValueType, BlockSideLength* BlockSideLength>
{
public:
    typedef typename swappable_block<ValueType, BlockSideLength* BlockSideLength>::internal_block_type internal_block_type;

    using swappable_block<ValueType, BlockSideLength* BlockSideLength>::get_internal_block;

    void fill_default()
    {
        // get_internal_block checks acquired
        internal_block_type& data = get_internal_block();
        #if STXXL_PARALLEL
        #pragma omp parallel for
        #endif
        for (int_type row = 0; row < int_type(BlockSideLength); ++row)
            for (int_type col = 0; col < int_type(BlockSideLength); ++col)
                data[row * BlockSideLength + col] = 0;
    }
};

//! External container for a (sub)matrix. Not intended for direct use.
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
//! \tparam BlockSideLength side length of a matrix block
//!
//! Stores blocks only, so all measures (height, width, row, col) are in blocks.
template <typename ValueType, unsigned BlockSideLength>
class swappable_block_matrix : public atomic_counted_object
{
public:
    typedef int_type size_type;
    typedef int_type elem_size_type;
    typedef block_scheduler<matrix_swappable_block<ValueType, BlockSideLength> > block_scheduler_type;
    typedef typename block_scheduler_type::swappable_block_identifier_type swappable_block_identifier_type;
    typedef std::vector<swappable_block_identifier_type> blocks_type;
    typedef matrix_local::matrix_operations<ValueType, BlockSideLength> Ops;

    block_scheduler_type& bs;

private:
    // assigning is not allowed
    swappable_block_matrix& operator = (const swappable_block_matrix& other);

protected:
    //! height of the matrix in blocks
    size_type height,
    //! width of the matrix in blocks
        width,
    //! height copied from supermatrix in blocks
        height_from_supermatrix,
    //! width copied from supermatrix in blocks
        width_from_supermatrix;
    //! the matrice's blocks in row-major
    blocks_type blocks;
    //! if the elements in each block are in col-major instead of row-major
    bool elements_in_blocks_transposed;

    //! get identifier of the block at (row, col)
    swappable_block_identifier_type & bl(const size_type row, const size_type col)
    { return blocks[row * width + col]; }

public:
    //! Create an empty swappable_block_matrix of given dimensions.
    swappable_block_matrix(block_scheduler_type& bs, const size_type height_in_blocks, const size_type width_in_blocks, const bool transposed = false)
        : bs(bs),
          height(height_in_blocks),
          width(width_in_blocks),
          height_from_supermatrix(0),
          width_from_supermatrix(0),
          blocks(height * width),
          elements_in_blocks_transposed(transposed)
    {
        for (size_type row = 0; row < height; ++row)
            for (size_type col = 0; col < width; ++col)
                bl(row, col) = bs.allocate_swappable_block();
    }

    //! Create swappable_block_matrix of given dimensions that
    //! represents the submatrix of supermatrix starting at (from_row_in_blocks, from_col_in_blocks).
    //!
    //! If supermatrix is not large enough, the submatrix is padded with empty blocks.
    //! The supermatrix must not be destructed or transposed before the submatrix is destructed.
    swappable_block_matrix(const swappable_block_matrix& supermatrix,
                           const size_type height_in_blocks, const size_type width_in_blocks,
                           const size_type from_row_in_blocks, const size_type from_col_in_blocks)
        : bs(supermatrix.bs),
          height(height_in_blocks),
          width(width_in_blocks),
          height_from_supermatrix(std::min(supermatrix.height - from_row_in_blocks, height)),
          width_from_supermatrix(std::min(supermatrix.width - from_col_in_blocks, width)),
          blocks(height * width),
          elements_in_blocks_transposed(supermatrix.elements_in_blocks_transposed)
    {
        for (size_type row = 0; row < height_from_supermatrix; ++row)
        {
            for (size_type col = 0; col < width_from_supermatrix; ++col)
                bl(row, col) = supermatrix.block(row + from_row_in_blocks, col + from_col_in_blocks);
            for (size_type col = width_from_supermatrix; col < width; ++col)
                bl(row, col) = bs.allocate_swappable_block();
        }
        for (size_type row = height_from_supermatrix; row < height; ++row)
            for (size_type col = 0; col < width; ++col)
                bl(row, col) = bs.allocate_swappable_block();
    }

    //! Create swappable_block_matrix that represents the combination matrix ul ur dl dr.
    //!
    //! The submatrices are assumed to be of fitting dimensions and equal transposition.
    //! The submatrices must not be destructed or transposed before the matrix is destructed.
    swappable_block_matrix(const swappable_block_matrix& ul, const swappable_block_matrix& ur,
                           const swappable_block_matrix& dl, const swappable_block_matrix& dr)
        : bs(ul.bs),
          height(ul.height + dl.height),
          width(ul.width + ur.width),
          height_from_supermatrix(height),
          width_from_supermatrix(width),
          blocks(height * width),
          elements_in_blocks_transposed(ul.elements_in_blocks_transposed)
    {
        for (size_type row = 0; row < ul.height; ++row)
        {
            for (size_type col = 0; col < ul.width; ++col)
                bl(row, col) = ul.block(row, col);
            for (size_type col = ul.width; col < width; ++col)
                bl(row, col) = ur.block(row, col - ul.width);
        }
        for (size_type row = ul.height; row < height; ++row)
        {
            for (size_type col = 0; col < ul.width; ++col)
                bl(row, col) = dl.block(row - ul.height, col);
            for (size_type col = ul.width; col < width; ++col)
                bl(row, col) = dr.block(row - ul.height, col - ul.width);
        }
    }

    swappable_block_matrix(const swappable_block_matrix& other)
        : atomic_counted_object(other),
          bs(other.bs),
          height(other.height),
          width(other.width),
          height_from_supermatrix(0),
          width_from_supermatrix(0),
          blocks(height * width),
          elements_in_blocks_transposed(false)
    {
        for (size_type row = 0; row < height; ++row)
            for (size_type col = 0; col < width; ++col)
                bl(row, col) = bs.allocate_swappable_block();
        // 0 + other is copying
        Ops::element_op(*this, other, typename Ops::addition());
    }

    ~swappable_block_matrix()
    {
        for (size_type row = 0; row < height_from_supermatrix; ++row)
        {
            for (size_type col = width_from_supermatrix; col < width; ++col)
                bs.free_swappable_block(bl(row, col));
        }
        for (size_type row = height_from_supermatrix; row < height; ++row)
            for (size_type col = 0; col < width; ++col)
                bs.free_swappable_block(bl(row, col));
    }

    static size_type block_index_from_elem(elem_size_type index)
    { return index / BlockSideLength; }

    static int_type elem_index_in_block_from_elem(elem_size_type index)
    { return index % BlockSideLength; }

    // regards transposed
    int_type elem_index_in_block_from_elem(elem_size_type row, elem_size_type col) const
    {
        return (is_transposed())
               ? row % BlockSideLength + col % BlockSideLength * BlockSideLength
               : row % BlockSideLength * BlockSideLength + col % BlockSideLength;
    }

    //! get identifier of the block at (row, col)
    const swappable_block_identifier_type & block(const size_type row, const size_type col) const
    { return blocks[row * width + col]; }

    //! get identifier of the block at (row, col)
    const swappable_block_identifier_type& operator () (const size_type row, const size_type col) const
    { return block(row, col); }

    const size_type & get_height() const
    { return height; }

    const size_type & get_width() const
    { return width; }

    //! if the elements inside the blocks are in transposed order i.e. column-major
    const bool & is_transposed() const
    { return elements_in_blocks_transposed; }

    void transpose()
    {
        // transpose matrix of blocks
        blocks_type bn(blocks.size());
        for (size_type row = 0; row < height; ++row)
            for (size_type col = 0; col < width; ++col)
                bn[col * height + row] = bl(row, col);
        bn.swap(blocks);
        // swap dimensions
        std::swap(height, width);
        std::swap(height_from_supermatrix, width_from_supermatrix);
        elements_in_blocks_transposed = ! elements_in_blocks_transposed;
    }

    void set_zero()
    {
        for (typename blocks_type::iterator it = blocks.begin(); it != blocks.end(); ++it)
            bs.deinitialize(*it);
    }
};

//! general iterator type that points to single elements inside a matrix
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
//! \tparam BlockSideLength side length of a matrix block
template <typename ValueType, unsigned BlockSideLength>
class matrix_iterator
{
protected:
    typedef matrix<ValueType, BlockSideLength> matrix_type;
    typedef typename matrix_type::swappable_block_matrix_type swappable_block_matrix_type;
    typedef typename matrix_type::block_scheduler_type block_scheduler_type;
    typedef typename block_scheduler_type::internal_block_type internal_block_type;
    typedef typename matrix_type::elem_size_type elem_size_type;
    typedef typename matrix_type::block_size_type block_size_type;

    template <typename VT, unsigned BSL>
    friend class matrix;

    template <typename VT, unsigned BSL>
    friend class const_matrix_iterator;

    matrix_type* m;
    elem_size_type current_row,          // \ both indices == -1 <=> empty iterator
        current_col;                     // /
    block_size_type current_block_row,
        current_block_col;
    internal_block_type* current_iblock; // NULL if block is not acquired

    void acquire_current_iblock()
    {
        if (! current_iblock)
            current_iblock = &m->data->bs.acquire(m->data->block(current_block_row, current_block_col));
    }

    void release_current_iblock()
    {
        if (current_iblock)
        {
            m->data->bs.release(m->data->block(current_block_row, current_block_col), true);
            current_iblock = 0;
        }
    }

    //! create iterator pointing to given row and col
    matrix_iterator(matrix_type& matrix, const elem_size_type start_row, const elem_size_type start_col)
        : m(&matrix),
          current_row(start_row),
          current_col(start_col),
          current_block_row(m->data->block_index_from_elem(start_row)),
          current_block_col(m->data->block_index_from_elem(start_col)),
          current_iblock(0) { }

    //! create empty iterator
    matrix_iterator(matrix_type& matrix)
        : m(&matrix),
          current_row(-1), // empty iterator
          current_col(-1),
          current_block_row(-1),
          current_block_col(-1),
          current_iblock(0) { }

    void set_empty()
    {
        release_current_iblock();
        current_row = -1;
        current_col = -1;
        current_block_row = -1;
        current_block_col = -1;
    }

public:
    matrix_iterator(const matrix_iterator& other)
        : m(other.m),
          current_row(other.current_row),
          current_col(other.current_col),
          current_block_row(other.current_block_row),
          current_block_col(other.current_block_col),
          current_iblock(0)
    {
        if (other.current_iblock)
            acquire_current_iblock();
    }

    matrix_iterator& operator = (const matrix_iterator& other)
    {
        set_pos(other.current_row, other.current_col);
        m = other.m;
        if (other.current_iblock)
            acquire_current_iblock();
        return *this;
    }

    ~matrix_iterator()
    { release_current_iblock(); }

    void set_row(const elem_size_type new_row)
    {
        const block_size_type new_block_row = m->data->block_index_from_elem(new_row);
        if (new_block_row != current_block_row)
        {
            release_current_iblock();
            current_block_row = new_block_row;
        }
        current_row = new_row;
    }

    void set_col(const elem_size_type new_col)
    {
        const block_size_type new_block_col = m->data->block_index_from_elem(new_col);
        if (new_block_col != current_block_col)
        {
            release_current_iblock();
            current_block_col = new_block_col;
        }
        current_col = new_col;
    }

    void set_pos(const elem_size_type new_row, const elem_size_type new_col)
    {
        const block_size_type new_block_row = m->data->block_index_from_elem(new_row),
            new_block_col = m->data->block_index_from_elem(new_col);
        if (new_block_col != current_block_col || new_block_row != current_block_row)
        {
            release_current_iblock();
            current_block_row = new_block_row;
            current_block_col = new_block_col;
        }
        current_row = new_row;
        current_col = new_col;
    }

    void set_pos(const std::pair<elem_size_type, elem_size_type> new_pos)
    { set_pos(new_pos.first, new_pos.second); }

    const elem_size_type & get_row() const
    { return current_row; }

    const elem_size_type & get_col() const
    { return current_col; }

    std::pair<elem_size_type, elem_size_type> get_pos() const
    { return std::make_pair(current_row, current_col); }

    bool empty() const
    { return current_row == -1 && current_col == -1; }

    operator bool () const
    { return ! empty(); }

    bool operator == (const matrix_iterator& other) const
    {
        return current_row == other.current_row && current_col == other.current_col && m == other.m;
    }

    //! Returns reference access to the element referenced by the iterator.
    //! The reference is only valid so long as the iterator is not moved.
    ValueType& operator * ()
    {
        acquire_current_iblock();
        return (*current_iblock)[m->data->elem_index_in_block_from_elem(current_row, current_col)];
    }
};

//! row-major iterator that points to single elements inside a matrix
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
//! \tparam BlockSideLength side length of a matrix block
template <typename ValueType, unsigned BlockSideLength>
class matrix_row_major_iterator : public matrix_iterator<ValueType, BlockSideLength>
{
protected:
    typedef matrix_iterator<ValueType, BlockSideLength> matrix_iterator_type;
    typedef typename matrix_iterator_type::matrix_type matrix_type;
    typedef typename matrix_iterator_type::elem_size_type elem_size_type;

    template <typename VT, unsigned BSL>
    friend class matrix;

    using matrix_iterator_type::m;
    using matrix_iterator_type::set_empty;

    //! create iterator pointing to given row and col
    matrix_row_major_iterator(matrix_type& matrix, const elem_size_type start_row, const elem_size_type start_col)
        : matrix_iterator_type(matrix, start_row, start_col) { }

    //! create empty iterator
    matrix_row_major_iterator(matrix_type& matrix)
        : matrix_iterator_type(matrix) { }

public:
    //! convert from matrix_iterator
    matrix_row_major_iterator(const matrix_iterator_type& matrix_iterator)
        : matrix_iterator_type(matrix_iterator) { }

    // Has to be not empty, else behavior is undefined.
    matrix_row_major_iterator& operator ++ ()
    {
        if (get_col() + 1 < m->get_width())
            // => not matrix_row_major_iterator the end of row, move right
            set_col(get_col() + 1);
        else if (get_row() + 1 < m->get_height())
            // => at end of row but not last row, move to beginning of next row
            set_pos(get_row() + 1, 0);
        else
            // => at end of matrix, set to empty-state
            set_empty();
        return *this;
    }

    // Has to be not empty, else behavior is undefined.
    matrix_row_major_iterator& operator -- ()
    {
        if (get_col() - 1 >= 0)
            // => not at the beginning of row, move left
            set_col(get_col() - 1);
        else if (get_row() - 1 >= 0)
            // => at beginning of row but not first row, move to end of previous row
            set_pos(get_row() - 1, m->get_width() - 1);
        else
            // => at beginning of matrix, set to empty-state
            set_empty();
        return *this;
    }

    using matrix_iterator_type::get_row;
    using matrix_iterator_type::get_col;
    using matrix_iterator_type::set_col;
    using matrix_iterator_type::set_pos;
};

//! column-major iterator that points to single elements inside a matrix
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
//! \tparam BlockSideLength side length of a matrix block
template <typename ValueType, unsigned BlockSideLength>
class matrix_col_major_iterator : public matrix_iterator<ValueType, BlockSideLength>
{
protected:
    typedef matrix_iterator<ValueType, BlockSideLength> matrix_iterator_type;
    typedef typename matrix_iterator_type::matrix_type matrix_type;
    typedef typename matrix_iterator_type::elem_size_type elem_size_type;

    template <typename VT, unsigned BSL>
    friend class matrix;

    using matrix_iterator_type::m;
    using matrix_iterator_type::set_empty;

    //! create iterator pointing to given row and col
    matrix_col_major_iterator(matrix_type& matrix, const elem_size_type start_row, const elem_size_type start_col)
        : matrix_iterator_type(matrix, start_row, start_col) { }

    //! create empty iterator
    matrix_col_major_iterator(matrix_type& matrix)
        : matrix_iterator_type(matrix) { }

public:
    //! convert from matrix_iterator
    matrix_col_major_iterator(const matrix_iterator_type& matrix_iterator)
        : matrix_iterator_type(matrix_iterator) { }

    // Has to be not empty, else behavior is undefined.
    matrix_col_major_iterator& operator ++ ()
    {
        if (get_row() + 1 < m->get_height())
            // => not at the end of col, move down
            set_row(get_row() + 1);
        else if (get_col() + 1 < m->get_width())
            // => at end of col but not last col, move to beginning of next col
            set_pos(0, get_col() + 1);
        else
            // => at end of matrix, set to empty-state
            set_empty();
        return *this;
    }

    // Has to be not empty, else behavior is undefined.
    matrix_col_major_iterator& operator -- ()
    {
        if (get_row() - 1 >= 0)
            // => not at the beginning of col, move up
            set_row(get_row() - 1);
        else if (get_col() - 1 >= 0)
            // => at beginning of col but not first col, move to end of previous col
            set_pos(m->get_height() - 1, get_col() - 1);
        else
            // => at beginning of matrix, set to empty-state
            set_empty();
        return *this;
    }

    using matrix_iterator_type::get_row;
    using matrix_iterator_type::get_col;
    using matrix_iterator_type::set_row;
    using matrix_iterator_type::set_pos;
};

//! general const_iterator type that points to single elements inside a matrix
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
//! \tparam BlockSideLength side length of a matrix block
template <typename ValueType, unsigned BlockSideLength>
class const_matrix_iterator
{
protected:
    typedef matrix<ValueType, BlockSideLength> matrix_type;
    typedef typename matrix_type::swappable_block_matrix_type swappable_block_matrix_type;
    typedef typename matrix_type::block_scheduler_type block_scheduler_type;
    typedef typename block_scheduler_type::internal_block_type internal_block_type;
    typedef typename matrix_type::elem_size_type elem_size_type;
    typedef typename matrix_type::block_size_type block_size_type;

    template <typename VT, unsigned BSL>
    friend class matrix;

    const matrix_type* m;
    elem_size_type current_row,          // \ both indices == -1 <=> empty iterator
        current_col;                     // /
    block_size_type current_block_row,
        current_block_col;
    internal_block_type* current_iblock; // NULL if block is not acquired

    void acquire_current_iblock()
    {
        if (! current_iblock)
            current_iblock = &m->data->bs.acquire(m->data->block(current_block_row, current_block_col));
    }

    void release_current_iblock()
    {
        if (current_iblock)
        {
            m->data->bs.release(m->data->block(current_block_row, current_block_col), false);
            current_iblock = 0;
        }
    }

    //! create iterator pointing to given row and col
    const_matrix_iterator(const matrix_type& matrix, const elem_size_type start_row, const elem_size_type start_col)
        : m(&matrix),
          current_row(start_row),
          current_col(start_col),
          current_block_row(m->data->block_index_from_elem(start_row)),
          current_block_col(m->data->block_index_from_elem(start_col)),
          current_iblock(0) { }

    //! create empty iterator
    const_matrix_iterator(const matrix_type& matrix)
        : m(&matrix),
          current_row(-1), // empty iterator
          current_col(-1),
          current_block_row(-1),
          current_block_col(-1),
          current_iblock(0) { }

    void set_empty()
    {
        release_current_iblock();
        current_row = -1;
        current_col = -1;
        current_block_row = -1;
        current_block_col = -1;
    }

public:
    const_matrix_iterator(const matrix_iterator<ValueType, BlockSideLength>& other)
        : m(other.m),
          current_row(other.current_row),
          current_col(other.current_col),
          current_block_row(other.current_block_row),
          current_block_col(other.current_block_col),
          current_iblock(0)
    {
        if (other.current_iblock)
            acquire_current_iblock();
    }

    const_matrix_iterator(const const_matrix_iterator& other)
        : m(other.m),
          current_row(other.current_row),
          current_col(other.current_col),
          current_block_row(other.current_block_row),
          current_block_col(other.current_block_col),
          current_iblock(0)
    {
        if (other.current_iblock)
            acquire_current_iblock();
    }

    const_matrix_iterator& operator = (const const_matrix_iterator& other)
    {
        set_pos(other.current_row, other.current_col);
        m = other.m;
        if (other.current_iblock)
            acquire_current_iblock();
        return *this;
    }

    ~const_matrix_iterator()
    { release_current_iblock(); }

    void set_row(const elem_size_type new_row)
    {
        const block_size_type new_block_row = m->data->block_index_from_elem(new_row);
        if (new_block_row != current_block_row)
        {
            release_current_iblock();
            current_block_row = new_block_row;
        }
        current_row = new_row;
    }

    void set_col(const elem_size_type new_col)
    {
        const block_size_type new_block_col = m->data->block_index_from_elem(new_col);
        if (new_block_col != current_block_col)
        {
            release_current_iblock();
            current_block_col = new_block_col;
        }
        current_col = new_col;
    }

    void set_pos(const elem_size_type new_row, const elem_size_type new_col)
    {
        const block_size_type new_block_row = m->data->block_index_from_elem(new_row),
            new_block_col = m->data->block_index_from_elem(new_col);
        if (new_block_col != current_block_col || new_block_row != current_block_row)
        {
            release_current_iblock();
            current_block_row = new_block_row;
            current_block_col = new_block_col;
        }
        current_row = new_row;
        current_col = new_col;
    }

    void set_pos(const std::pair<elem_size_type, elem_size_type> new_pos)
    { set_pos(new_pos.first, new_pos.second); }

    const elem_size_type & get_row() const
    { return current_row; }

    const elem_size_type & get_col() const
    { return current_col; }

    std::pair<elem_size_type, elem_size_type> get_pos() const
    { return std::make_pair(current_row, current_col); }

    bool empty() const
    { return current_row == -1 && current_col == -1; }

    operator bool () const
    { return ! empty(); }

    bool operator == (const const_matrix_iterator& other) const
    {
        return current_row == other.current_row && current_col == other.current_col && m == other.m;
    }

    //! Returns reference access to the element referenced by the iterator.
    //! The reference is only valid so long as the iterator is not moved.
    const ValueType& operator * ()
    {
        acquire_current_iblock();
        return (*current_iblock)[m->data->elem_index_in_block_from_elem(current_row, current_col)];
    }
};

//! row-major const_iterator that points to single elements inside a matrix
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
//! \tparam BlockSideLength side length of a matrix block
template <typename ValueType, unsigned BlockSideLength>
class const_matrix_row_major_iterator : public const_matrix_iterator<ValueType, BlockSideLength>
{
protected:
    typedef const_matrix_iterator<ValueType, BlockSideLength> const_matrix_iterator_type;
    typedef typename const_matrix_iterator_type::matrix_type matrix_type;
    typedef typename const_matrix_iterator_type::elem_size_type elem_size_type;

    template <typename VT, unsigned BSL>
    friend class matrix;

    using const_matrix_iterator_type::m;
    using const_matrix_iterator_type::set_empty;

    //! create iterator pointing to given row and col
    const_matrix_row_major_iterator(const matrix_type& matrix, const elem_size_type start_row, const elem_size_type start_col)
        : const_matrix_iterator_type(matrix, start_row, start_col) { }

    //! create empty iterator
    const_matrix_row_major_iterator(const matrix_type& matrix)
        : const_matrix_iterator_type(matrix) { }

public:
    //! convert from matrix_iterator
    const_matrix_row_major_iterator(const const_matrix_row_major_iterator& matrix_iterator)
        : const_matrix_iterator_type(matrix_iterator) { }

    //! convert from matrix_iterator
    const_matrix_row_major_iterator(const const_matrix_iterator_type& matrix_iterator)
        : const_matrix_iterator_type(matrix_iterator) { }

    // Has to be not empty, else behavior is undefined.
    const_matrix_row_major_iterator& operator ++ ()
    {
        if (get_col() + 1 < m->get_width())
            // => not matrix_row_major_iterator the end of row, move right
            set_col(get_col() + 1);
        else if (get_row() + 1 < m->get_height())
            // => at end of row but not last row, move to beginning of next row
            set_pos(get_row() + 1, 0);
        else
            // => at end of matrix, set to empty-state
            set_empty();
        return *this;
    }

    // Has to be not empty, else behavior is undefined.
    const_matrix_row_major_iterator& operator -- ()
    {
        if (get_col() - 1 >= 0)
            // => not at the beginning of row, move left
            set_col(get_col() - 1);
        else if (get_row() - 1 >= 0)
            // => at beginning of row but not first row, move to end of previous row
            set_pos(get_row() - 1, m->get_width() - 1);
        else
            // => at beginning of matrix, set to empty-state
            set_empty();
        return *this;
    }

    using const_matrix_iterator_type::get_row;
    using const_matrix_iterator_type::get_col;
    using const_matrix_iterator_type::set_col;
    using const_matrix_iterator_type::set_pos;
};

//! column-major const_iterator that points to single elements inside a matrix
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
//! \tparam BlockSideLength side length of a matrix block
template <typename ValueType, unsigned BlockSideLength>
class const_matrix_col_major_iterator : public const_matrix_iterator<ValueType, BlockSideLength>
{
protected:
    typedef const_matrix_iterator<ValueType, BlockSideLength> const_matrix_iterator_type;
    typedef typename const_matrix_iterator_type::matrix_type matrix_type;
    typedef typename const_matrix_iterator_type::elem_size_type elem_size_type;

    template <typename VT, unsigned BSL>
    friend class matrix;

    using const_matrix_iterator_type::m;
    using const_matrix_iterator_type::set_empty;

    //! create iterator pointing to given row and col
    const_matrix_col_major_iterator(const matrix_type& matrix, const elem_size_type start_row, const elem_size_type start_col)
        : const_matrix_iterator_type(matrix, start_row, start_col) { }

    //! create empty iterator
    const_matrix_col_major_iterator(const matrix_type& matrix)
        : const_matrix_iterator_type(matrix) { }

public:
    //! convert from matrix_iterator
    const_matrix_col_major_iterator(const matrix_iterator<ValueType, BlockSideLength>& matrix_iterator)
        : const_matrix_iterator_type(matrix_iterator) { }

    //! convert from matrix_iterator
    const_matrix_col_major_iterator(const const_matrix_iterator_type& matrix_iterator)
        : const_matrix_iterator_type(matrix_iterator) { }

    // Has to be not empty, else behavior is undefined.
    const_matrix_col_major_iterator& operator ++ ()
    {
        if (get_row() + 1 < m->get_height())
            // => not at the end of col, move down
            set_row(get_row() + 1);
        else if (get_col() + 1 < m->get_width())
            // => at end of col but not last col, move to beginning of next col
            set_pos(0, get_col() + 1);
        else
            // => at end of matrix, set to empty-state
            set_empty();
        return *this;
    }

    // Has to be not empty, else behavior is undefined.
    const_matrix_col_major_iterator& operator -- ()
    {
        if (get_row() - 1 >= 0)
            // => not at the beginning of col, move up
            set_row(get_row() - 1);
        else if (get_col() - 1 >= 0)
            // => at beginning of col but not first col, move to end of previous col
            set_pos(m->get_height() - 1, get_col() - 1);
        else
            // => at beginning of matrix, set to empty-state
            set_empty();
        return *this;
    }

    using const_matrix_iterator_type::get_row;
    using const_matrix_iterator_type::get_col;
    using const_matrix_iterator_type::set_row;
    using const_matrix_iterator_type::set_pos;
};

//! External matrix container. \n
//! <b> Introduction </b> to matrix container: see \ref tutorial_matrix tutorial. \n
//! <b> Design and Internals </b> of matrix container: see \ref design_matrix.
//!
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
//! \tparam BlockSideLength side length of a matrix block
//!
//! Divides the matrix in square submatrices (blocks).
//! Blocks can be swapped individually to and from external memory.
//! They are only swapped if necessary to minimize I/O.
template <typename ValueType, unsigned BlockSideLength>
class matrix
{
protected:
    typedef matrix<ValueType, BlockSideLength> matrix_type;
    typedef swappable_block_matrix<ValueType, BlockSideLength> swappable_block_matrix_type;
    typedef counting_ptr<swappable_block_matrix_type> swappable_block_matrix_pointer_type;
    typedef typename swappable_block_matrix_type::block_scheduler_type block_scheduler_type;
    typedef typename swappable_block_matrix_type::size_type block_size_type;
    typedef typename swappable_block_matrix_type::elem_size_type elem_size_type;
    typedef matrix_local::matrix_operations<ValueType, BlockSideLength> Ops;
    typedef matrix_swappable_block<ValueType, BlockSideLength> swappable_block_type;

public:
    typedef matrix_iterator<ValueType, BlockSideLength> iterator;
    typedef const_matrix_iterator<ValueType, BlockSideLength> const_iterator;
    typedef matrix_row_major_iterator<ValueType, BlockSideLength> row_major_iterator;
    typedef matrix_col_major_iterator<ValueType, BlockSideLength> col_major_iterator;
    typedef const_matrix_row_major_iterator<ValueType, BlockSideLength> const_row_major_iterator;
    typedef const_matrix_col_major_iterator<ValueType, BlockSideLength> const_col_major_iterator;
    typedef column_vector<ValueType> column_vector_type;
    typedef row_vector<ValueType> row_vector_type;

protected:
    template <typename VT, unsigned BSL>
    friend class matrix_iterator;

    template <typename VT, unsigned BSL>
    friend class const_matrix_iterator;

    elem_size_type height, width;

    swappable_block_matrix_pointer_type data;

public:
    //! \name Constructors/Destructors
    //! \{

    //! Creates a new matrix of given dimensions. Elements' values are set to zero.
    //! \param bs block scheduler used
    //! \param height height of the created matrix
    //! \param width width of the created matrix
    matrix(block_scheduler_type& bs,
           const elem_size_type height, const elem_size_type width)
        : height(height),
          width(width),
          data(
              new swappable_block_matrix_type(
                  bs,
                  div_ceil(height, BlockSideLength),
                  div_ceil(width, BlockSideLength))
              )
    { }

    matrix(block_scheduler_type& bs,
           const column_vector_type& left, const row_vector_type& right)
        : height((elem_size_type)left.size()),
          width((elem_size_type)right.size()),
          data(
              new swappable_block_matrix_type(
                  bs,
                  div_ceil(height, BlockSideLength),
                  div_ceil(width, BlockSideLength))
              )
    { Ops::recursive_matrix_from_vectors(*data, left, right); }

    ~matrix() { }
    //! \}

    //! \name Capacity
    //! \{
    const elem_size_type & get_height() const
    { return height; }

    const elem_size_type & get_width() const
    { return width; }
    //! \}

    //! \name Iterators
    //! \{
    iterator begin()
    {
        data.unify();
        return iterator(*this, 0, 0);
    }
    const_iterator begin() const
    { return const_iterator(*this, 0, 0); }
    const_iterator cbegin() const
    { return const_iterator(*this, 0, 0); }
    iterator end()
    {
        data.unify();
        return iterator(*this);
    }
    const_iterator end() const
    { return const_iterator(*this); }
    const_iterator cend() const
    { return const_iterator(*this); }
    const_iterator operator () (const elem_size_type row, const elem_size_type col) const
    { return const_iterator(*this, row, col); }

    iterator operator () (const elem_size_type row, const elem_size_type col)
    {
        data.unify();
        return iterator(*this, row, col);
    }
    //! \}

    //! \name Modifiers
    //! \{
    void transpose()
    {
        data.unify();
        data->transpose();
        std::swap(height, width);
    }

    void set_zero()
    {
        if (data.unique())
            data->set_zero();
        else
            data = new swappable_block_matrix_type
                       (data->bs, div_ceil(height, BlockSideLength), div_ceil(width, BlockSideLength));
    }
    //! \}

    //! \name Operations
    //! \{
    matrix_type operator + (const matrix_type& right) const
    {
        assert(height == right.height && width == right.width);
        matrix_type res(data->bs, height, width);
        Ops::element_op(*res.data, *data, *right.data, typename Ops::addition()); // more efficient than copying this and then adding right
        return res;
    }

    matrix_type operator - (const matrix_type& right) const
    {
        assert(height == right.height && width == right.width);
        matrix_type res(data->bs, height, width);
        Ops::element_op(*res.data, *data, *right.data, typename Ops::subtraction()); // more efficient than copying this and then subtracting right
        return res;
    }

    matrix_type operator * (const matrix_type& right) const
    { return multiply(right); }

    matrix_type operator * (const ValueType scalar) const
    {
        matrix_type res(data->bs, height, width);
        Ops::element_op(*res.data, *data, typename Ops::scalar_multiplication(scalar));
        return res;
    }

    matrix_type& operator += (const matrix_type& right)
    {
        assert(height == right.height && width == right.width);
        data.unify();
        Ops::element_op(*data, *right.data, typename Ops::addition());
        return *this;
    }

    matrix_type& operator -= (const matrix_type& right)
    {
        assert(height == right.height && width == right.width);
        data.unify();
        Ops::element_op(*data, *right.data, typename Ops::subtraction());
        return *this;
    }

    matrix_type& operator *= (const matrix_type& right)
    { return *this = operator * (right); } // implicitly unifies by constructing a result-matrix

    matrix_type& operator *= (const ValueType scalar)
    {
        data.unify();
        Ops::element_op(*data, typename Ops::scalar_multiplication(scalar));
        return *this;
    }

    column_vector_type operator * (const column_vector_type& right) const
    {
        assert(elem_size_type(right.size()) == width);
        column_vector_type res(height);
        res.set_zero();
        Ops::recursive_matrix_col_vector_multiply_and_add(*data, right, res);
        return res;
    }

    row_vector_type multiply_from_left(const row_vector_type& left) const
    {
        assert(elem_size_type(left.size()) == height);
        row_vector_type res(width);
        res.set_zero();
        Ops::recursive_matrix_row_vector_multiply_and_add(left, *data, res);
        return res;
    }

    //! multiply with another matrix
    //! \param right matrix to multiply with
    //! \param multiplication_algorithm allows to choose the applied algorithm
    //! \param scheduling_algorithm  allows to choose the applied algorithm
    //!
    //! Available algorithms are: \n
    //!    0: naive_multiply_and_add (I/O inefficient, slow) \n
    //!    1: recursive_multiply_and_add (recommended, default, stable time and I/O complexity) \n
    //!    2: strassen_winograd_multiply_and_add (sometimes fast but unstable time and I/O complexity) \n
    //!    3: multi_level_strassen_winograd_multiply_and_add (sometimes fast but unstable time and I/O complexity) \n
    //!    4: strassen_winograd_multiply, optimized pre- and postadditions (sometimes fast but unstable time and I/O complexity) \n
    //!    5: strassen_winograd_multiply_and_add_interleaved, optimized preadditions (sometimes fast but unstable time and I/O complexity) \n
    //!    6: multi_level_strassen_winograd_multiply_and_add_block_grained (sometimes fast but unstable time and I/O complexity)
    matrix_type multiply(const matrix_type& right, const int_type multiplication_algorithm = 1, const int_type scheduling_algorithm = 2) const
    {
        assert(width == right.height);
        assert(&data->bs == &right.data->bs);
        matrix_type res(data->bs, height, right.width);

        if (scheduling_algorithm > 0)
        {
            // all offline algos need a simulation-run
            delete data->bs.switch_algorithm_to(
                new block_scheduler_algorithm_simulation<swappable_block_type>(data->bs)
                );
            switch (multiplication_algorithm)
            {
            case 0:
                Ops::naive_multiply_and_add(*data, *right.data, *res.data);
                break;
            case 1:
                Ops::recursive_multiply_and_add(*data, *right.data, *res.data);
                break;
            case 2:
                Ops::strassen_winograd_multiply_and_add(*data, *right.data, *res.data);
                break;
            case 3:
                Ops::multi_level_strassen_winograd_multiply_and_add(*data, *right.data, *res.data);
                break;
            case 4:
                Ops::strassen_winograd_multiply(*data, *right.data, *res.data);
                break;
            case 5:
                Ops::strassen_winograd_multiply_and_add_interleaved(*data, *right.data, *res.data);
                break;
            case 6:
                Ops::multi_level_strassen_winograd_multiply_and_add_block_grained(*data, *right.data, *res.data);
                break;
            default:
                STXXL_ERRMSG("invalid multiplication-algorithm number");
                break;
            }
        }
        switch (scheduling_algorithm)
        {
        case 0:
            delete data->bs.switch_algorithm_to(
                new block_scheduler_algorithm_online_lru<swappable_block_type>(data->bs)
                );
            break;
        case 1:
            delete data->bs.switch_algorithm_to(
                new block_scheduler_algorithm_offline_lfd<swappable_block_type>(data->bs)
                );
            break;
        case 2:
            delete data->bs.switch_algorithm_to(
                new block_scheduler_algorithm_offline_lru_prefetching<swappable_block_type>(data->bs)
                );
            break;
        default:
            STXXL_ERRMSG("invalid scheduling-algorithm number");
        }
        switch (multiplication_algorithm)
        {
        case 0:
            Ops::naive_multiply_and_add(*data, *right.data, *res.data);
            break;
        case 1:
            Ops::recursive_multiply_and_add(*data, *right.data, *res.data);
            break;
        case 2:
            Ops::strassen_winograd_multiply_and_add(*data, *right.data, *res.data);
            break;
        case 3:
            Ops::multi_level_strassen_winograd_multiply_and_add(*data, *right.data, *res.data);
            break;
        case 4:
            Ops::strassen_winograd_multiply(*data, *right.data, *res.data);
            break;
        case 5:
            Ops::strassen_winograd_multiply_and_add_interleaved(*data, *right.data, *res.data);
            break;
        case 6:
            Ops::multi_level_strassen_winograd_multiply_and_add_block_grained(*data, *right.data, *res.data);
            break;
        default:
            STXXL_ERRMSG("invalid multiplication-algorithm number");
            break;
        }
        delete data->bs.switch_algorithm_to(
            new block_scheduler_algorithm_online_lru<swappable_block_type>(data->bs)
            );
        return res;
    }

    //! Use internal memory multiplication. Designated for testing. May exceed memory limitations.
    matrix_type multiply_internal(const matrix_type& right, const int_type scheduling_algorithm = 2) const
    {
        assert(width == right.height);
        assert(&data->bs == &right.data->bs);
        matrix_type res(data->bs, height, right.width);

        if (scheduling_algorithm > 0)
        {
            // all offline algos need a simulation-run
            delete data->bs.switch_algorithm_to(
                new block_scheduler_algorithm_simulation<swappable_block_type>(data->bs)
                );
            multiply_internal(right, res);
        }
        switch (scheduling_algorithm)
        {
        case 0:
            delete data->bs.switch_algorithm_to(
                new block_scheduler_algorithm_online_lru<swappable_block_type>(data->bs)
                );
            break;
        case 1:
            delete data->bs.switch_algorithm_to(
                new block_scheduler_algorithm_offline_lfd<swappable_block_type>(data->bs)
                );
            break;
        case 2:
            delete data->bs.switch_algorithm_to(
                new block_scheduler_algorithm_offline_lru_prefetching<swappable_block_type>(data->bs)
                );
            break;
        default:
            STXXL_ERRMSG("invalid scheduling-algorithm number");
        }
        multiply_internal(right, res);
        delete data->bs.switch_algorithm_to(
            new block_scheduler_algorithm_online_lru<swappable_block_type>(data->bs)
            );
        return res;
    }
    //! \}

protected:
    void multiply_internal(const matrix_type& right, matrix_type& res) const
    {
        ValueType* A = new ValueType[height * width];
        ValueType* B = new ValueType[right.height * right.width];
        ValueType* C = new ValueType[res.height * res.width];
        ValueType* vit;
        vit = A;
        for (const_row_major_iterator mit = cbegin(); mit != cend(); ++mit, ++vit)
            *vit = *mit;
        vit = B;
        for (const_row_major_iterator mit = right.cbegin(); mit != right.cend(); ++mit, ++vit)
            *vit = *mit;
        if (! res.data->bs.is_simulating())
        {
#if STXXL_BLAS
            gemm_wrapper(height, width, res.width,
                         ValueType(1), false, A,
                         false, B,
                         ValueType(0), false, C);
#else
            assert(false /* internal multiplication is only available for testing with blas */);
#endif
        }
        vit = C;
        for (row_major_iterator mit = res.begin(); mit != res.end(); ++mit, ++vit)
            *mit = *vit;
        delete[] A;
        delete[] B;
        delete[] C;
    }
};

//! \}

STXXL_END_NAMESPACE

#endif // !STXXL_CONTAINERS_MATRIX_HEADER
// vim: et:ts=4:sw=4