matrix_recursator.hpp

矩阵运算源码最新版本

	: my_sub_matrix(constructor_helper(matrix)), my_bound(outer_bound(matrix)),
	  splitter(my_sub_matrix)
    {
      if (bound == 0)
	my_bound= outer_bound(matrix);
      else {
	assert(is_power_of_2(bound));
	assert(bound >= matrix.num_rows() && bound >= matrix.num_cols());
	my_bound= bound;
      }
    }

    // Sub-matrices are copied directly
    // explicit matrix_recursator(sub_matrix_type sub_matrix) : my_sub_matrix(sub_matrix) {}
    
    sub_matrix_type& get_value()
    {
	return my_sub_matrix;
    }

    sub_matrix_type const& get_value() const
    {
	return my_sub_matrix;
    }

    // Returning quadrants for non-const recursator

    self north_west()
    {
	sub_matrix_type sm(sub_matrix(my_sub_matrix, my_sub_matrix.begin_row(), splitter.row_split(),
				      my_sub_matrix.begin_col(), splitter.col_split()));
	self tmp(sm, my_bound / 2);
	return tmp;
    }

    self south_west()
    {
	sub_matrix_type sm(sub_matrix(my_sub_matrix, splitter.row_split(), my_sub_matrix.end_row(), 
				      my_sub_matrix.begin_col(), splitter.col_split()));
	self tmp(sm, my_bound / 2);
	return tmp;
    }

    self north_east()
    {
	sub_matrix_type sm(sub_matrix(my_sub_matrix, my_sub_matrix.begin_row(), splitter.row_split(),
				      splitter.col_split(), my_sub_matrix.end_col()));
	self tmp(sm, my_bound / 2);
	return tmp;
    }

    self south_east()
    {
	sub_matrix_type sm(sub_matrix(my_sub_matrix, splitter.row_split(), my_sub_matrix.end_row(), 
				      splitter.col_split(), my_sub_matrix.end_col()));
	self tmp(sm, my_bound / 2);
	return tmp;
    }

    // Returning quadrants for const recursator

    self const north_west() const
    {
	sub_matrix_type sm(sub_matrix(const_cast<self*>(this)->my_sub_matrix, my_sub_matrix.begin_row(), splitter.row_split(),
				      my_sub_matrix.begin_col(), splitter.col_split()));
	self tmp(sm, my_bound / 2);
	return tmp;
    }

    self const south_west() const 
    {
	sub_matrix_type sm(sub_matrix(const_cast<self*>(this)->my_sub_matrix, splitter.row_split(), my_sub_matrix.end_row(), 
				      my_sub_matrix.begin_col(), splitter.col_split()));
	self tmp(sm, my_bound / 2);
	return tmp;
    }

    self const north_east() const 
    {
	sub_matrix_type sm(sub_matrix(const_cast<self*>(this)->my_sub_matrix, my_sub_matrix.begin_row(), splitter.row_split(),
				      splitter.col_split(), my_sub_matrix.end_col()));
	self tmp(sm, my_bound / 2);
	return tmp;
    }

    self const south_east() const 
    {
	sub_matrix_type sm(sub_matrix(const_cast<self*>(this)->my_sub_matrix, splitter.row_split(), my_sub_matrix.end_row(), 
				      splitter.col_split(), my_sub_matrix.end_col()));
	self tmp(sm, my_bound / 2);
	return tmp;
    }

    // Checking whether a quadrant is empty

    // For completeness
    bool north_west_empty() const
    {
	return false;
    }

    bool north_east_empty() const
    {
	return splitter.col_split() == my_sub_matrix.end_col();
    }

    bool south_west_empty() const
    {
	return splitter.row_split() == my_sub_matrix.end_row();
    }

    bool south_east_empty() const
    {
	return splitter.row_split() == my_sub_matrix.end_row() 
	       || splitter.col_split() == my_sub_matrix.end_col();
    }

    bool is_empty() const
    {
	return my_sub_matrix.begin_row() == my_sub_matrix.end_row()
	       || my_sub_matrix.begin_col() == my_sub_matrix.end_col();
    }

#if 0
    bool is_leaf() const
    {
	return my_sub_matrix.num_rows() < 2 || my_sub_matrix.num_cols() < 2;
    }
#endif

    size_type bound() const
    {
	assert(my_bound >= my_sub_matrix.num_rows() && my_bound >= my_sub_matrix.num_cols());
	return my_bound;
    }

    template <typename R1, typename R2> friend void equalize_depth (R1&, R2&);   
    template <typename R1, typename R2, typename R3> friend void equalize_depth (R1&, R2&, R3&);

  protected:
    sub_matrix_type     my_sub_matrix;
    size_type           my_bound;
    splitter_type       splitter;
};


template <typename Recursator1, typename Recursator2>
void inline equalize_depth(Recursator1& r1, Recursator2& r2)
{
    typename Recursator1::size_type max_bound= std::max(r1.bound(), r2.bound());
    r1.my_bound= max_bound;
    r2.my_bound= max_bound;
}

template <typename Recursator1, typename Recursator2, typename Recursator3>
void inline equalize_depth(Recursator1& r1, Recursator2& r2, Recursator3& r3)
{
    typename Recursator1::size_type max_bound= std::max(std::max(r1.bound(), r2.bound()), r3.bound());
    r1.my_bound= max_bound;
    r2.my_bound= max_bound;
    r3.my_bound= max_bound;
}


} // namespace recursion

/*! Import matrix_recursator into the mtl namespace
    \sa recursion::matrix_recursator, \ref rec_intro "recursion intro"
**/
using recursion::matrix_recursator;

// Define free functions (from member functions)

/*! Compute the north-west quadrant of a recursator (i.e. its referred matrix).
    The result is itself a recursator.
    \sa \ref rec_intro "recursion intro"
**/
template <typename Matrix>
matrix_recursator<Matrix> inline north_west(const matrix_recursator<Matrix>& rec)
{
    return rec.north_west();
}

/*! Compute the north-east quadrant of a recursator (i.e. its referred matrix).
    The result is itself a recursator.
    \sa \ref rec_intro "recursion intro"
**/
template <typename Matrix>
matrix_recursator<Matrix> inline north_east(const matrix_recursator<Matrix>& rec)
{
    return rec.north_east();
}

/*! Compute the south-west quadrant of a recursator (i.e. its referred matrix).
    The result is itself a recursator.
    \sa \ref rec_intro "recursion intro"
**/
template <typename Matrix>
matrix_recursator<Matrix> inline south_west(const matrix_recursator<Matrix>& rec)
{
    return rec.south_west();
}

/*! Compute the south-east quadrant of a recursator (i.e. its referred matrix).
    The result is itself a recursator.
    \sa \ref rec_intro "recursion intro"
**/
template <typename Matrix>
matrix_recursator<Matrix> inline south_east(const matrix_recursator<Matrix>& rec)
{
    return rec.south_east();
}


/*! Check if a recursator (i.e. its referred matrix) is empty.
    \sa \ref rec_intro "recursion intro"
**/
template <typename Matrix>
bool inline is_empty(const matrix_recursator<Matrix>& rec)
{
    return rec.is_empty();
}

/*! Check if a recursator (i.e. its referred matrix) fills the
    entire block, i.e. if the number of rows and columns are both
    equal to the virtual bound.
    \sa \ref rec_intro "recursion intro"
**/
template <typename Matrix>
bool inline is_full(const matrix_recursator<Matrix>& rec)
{
    int nr= rec.num_rows(), nc= rec.num_cols(), b= rec.bound();
    return rec.num_rows() == rec.bound() && rec.num_cols() == rec.bound();
}

/*! The number of rows that a sub-matrix would have if it was constructed.
    \sa \ref rec_intro "recursion intro"
**/
template <typename Matrix>
typename matrix_recursator<Matrix>::size_type
inline num_rows(const matrix_recursator<Matrix>& rec)
{
    return rec.num_rows();
}

/*! The number of columns that a sub-matrix would have if it was constructed.
    \sa \ref rec_intro "recursion intro"
**/
template <typename Matrix>
typename matrix_recursator<Matrix>::size_type
inline num_cols(const matrix_recursator<Matrix>& rec)
{
    return rec.num_cols();
}

/*! The number of elements (rows times columns) that a sub-matrix would have if it was constructed.
    \sa \ref rec_intro "recursion intro"
**/
template <typename Matrix>
typename matrix_recursator<Matrix>::size_type
inline size(const matrix_recursator<Matrix>& rec)
{
    return num_rows(rec) * num_cols(rec);
}


} // namespace mtl

#endif // MTL_MATRIX_RECURATOR_INCLUDE