/usr/include/stxxl/bits/containers
Edit: /usr/include/stxxl/bits/containers/matrix.h (51840B)
/***************************************************************************
* include/stxxl/bits/containers/matrix.h
*
* Part of the STXXL. See http://stxxl.sourceforge.net
*
* Copyright (C) 2010-2011 Raoul Steffen
*
* 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
#include
#include
#include
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
class matrix;
//! external column-vector container for matrix multiplication
//! \tparam ValueType type of contained objects (POD with no references to internal memory)
template
class column_vector : public vector
{
public:
typedef vector 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
class row_vector : public vector
{
public:
typedef vector 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
row_vector operator * (const matrix& right) const
{ return right.multiply_from_left(*this); }
ValueType operator * (const column_vector& 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
class matrix_swappable_block : public swappable_block
{
public:
typedef typename swappable_block::internal_block_type internal_block_type;
using swappable_block::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
class swappable_block_matrix : public atomic_counted_object
{
public:
typedef int_type size_type;
typedef int_type elem_size_type;
typedef block_scheduler > block_scheduler_type;
typedef typename block_scheduler_type::swappable_block_identifier_type swappable_block_identifier_type;
typedef std::vector blocks_type;
typedef matrix_local::matrix_operations 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
class matrix_iterator
{
protected:
typedef matrix 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
friend class matrix;
template
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 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 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
class matrix_row_major_iterator : public matrix_iterator
{
protected:
typedef matrix_iterator matrix_iterator_type;
typedef typename matrix_iterator_type::matrix_type matrix_type;
typedef typename matrix_iterator_type::elem_size_type elem_size_type;
template
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
class matrix_col_major_iterator : public matrix_iterator
{
protected:
typedef matrix_iterator matrix_iterator_type;
typedef typename matrix_iterator_type::matrix_type matrix_type;
typedef typename matrix_iterator_type::elem_size_type elem_size_type;
template
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
class const_matrix_iterator
{
protected:
typedef matrix 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
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& 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 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 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
class const_matrix_row_major_iterator : public const_matrix_iterator
{
protected:
typedef const_matrix_iterator 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
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
class const_matrix_col_major_iterator : public const_matrix_iterator
{
protected:
typedef const_matrix_iterator 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
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& 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
//! Introduction to matrix container: see \ref tutorial_matrix tutorial. \n
//! Design and Internals 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
class matrix
{
protected:
typedef matrix matrix_type;
typedef swappable_block_matrix swappable_block_matrix_type;
typedef counting_ptr 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 Ops;
typedef matrix_swappable_block swappable_block_type;
public:
typedef matrix_iterator iterator;
typedef const_matrix_iterator const_iterator;
typedef matrix_row_major_iterator row_major_iterator;
typedef matrix_col_major_iterator col_major_iterator;
typedef const_matrix_row_major_iterator const_row_major_iterator;
typedef const_matrix_col_major_iterator const_col_major_iterator;
typedef column_vector column_vector_type;
typedef row_vector row_vector_type;
protected:
template
friend class matrix_iterator;
template
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(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(data->bs)
);
break;
case 1:
delete data->bs.switch_algorithm_to(
new block_scheduler_algorithm_offline_lfd(data->bs)
);
break;
case 2:
delete data->bs.switch_algorithm_to(
new block_scheduler_algorithm_offline_lru_prefetching(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(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(data->bs)
);
multiply_internal(right, res);
}
switch (scheduling_algorithm)
{
case 0:
delete data->bs.switch_algorithm_to(
new block_scheduler_algorithm_online_lru(data->bs)
);
break;
case 1:
delete data->bs.switch_algorithm_to(
new block_scheduler_algorithm_offline_lfd(data->bs)
);
break;
case 2:
delete data->bs.switch_algorithm_to(
new block_scheduler_algorithm_offline_lru_prefetching(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(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