matrix_recursator.hpp

矩阵运算源码最新版本

// Software License for MTL
// 
// Copyright (c) 2007 The Trustees of Indiana University. All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
// 
// This file is part of the Matrix Template Library
// 
// See also license.mtl.txt in the distribution.

#ifndef MTL_MATRIX_RECURATOR_INCLUDE
#define MTL_MATRIX_RECURATOR_INCLUDE

#include <cmath>
#include <boost/numeric/mtl/recursion/utility.hpp>
#include <boost/numeric/mtl/operation/sub_matrix.hpp>
#include <boost/numeric/mtl/matrix/transposed_view.hpp>
#include <boost/numeric/mtl/recursion/dim_splitter.hpp>
#include <boost/numeric/mtl/operation/print_matrix.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>


namespace mtl { namespace recursion {


template <typename Recursator1, typename Recursator2>
void inline equalize_depth(Recursator1& r1, Recursator2& r2);

template <typename Recursator1, typename Recursator2, typename Recursator3>
void inline equalize_depth(Recursator1& r1, Recursator2& r2, Recursator3& r3);


/*! Class for matrix recursator

    How to use this class is described in the \ref rec_intro "recursion introduction".

    \sa \ref mtl::north_west, \ref mtl::north_east, 
    \ref mtl::south_west, \ref mtl::south_east, 
    \ref mtl::is_empty(const matrix_recursator<Matrix>&), 
    \ref mtl::is_full(const matrix_recursator<Matrix>&), 
    \ref mtl::num_rows(const matrix_recursator<Matrix>&), 
    \ref mtl::num_cols(const matrix_recursator<Matrix>&),
    \ref mtl::size(const matrix_recursator<Matrix>&)
**/
template <typename Matrix>
struct matrix_recursator
{
    typedef matrix_recursator                                     self;
    typedef Matrix                                                matrix_type;
    typedef typename sub_matrix_t<Matrix>::sub_matrix_type        sub_matrix_type;
    typedef typename sub_matrix_t<Matrix>::const_sub_matrix_type  const_sub_matrix_type;
    typedef typename Matrix::size_type                            size_type;
    typedef outer_bound_splitter<self>                            splitter_type;

private:
    
    template <typename MatrixType> // why was it templated ???
    sub_matrix_type constructor_helper(MatrixType const& matrix)
    {
	return sub_matrix(matrix, matrix.begin_row(), matrix.end_row(),
			  matrix.begin_col(), matrix.end_col());
    }

    // For views without own data, we need to generate a new sub_matrix as shared_ptr
    template <typename MatrixType>
    sub_matrix_type constructor_helper(transposed_view<MatrixType> const& view)
    {
	typedef typename sub_matrix_t<MatrixType>::sub_matrix_type   ref_sub_type;
	typedef boost::shared_ptr<ref_sub_type>                  pointer_type;
	typedef typename transposed_view<MatrixType>::other      ref_type;

	// Submatrix of referred matrix, colums and rows interchanged
	// Create a submatrix, whos address will be kept by transposed_view
	pointer_type p(new ref_sub_type(sub_matrix(const_cast<ref_type&>(view.ref), view.begin_col(), view.end_col(), 
						   view.begin_row(), view.end_row())));
	return sub_matrix_type(p); 
    }

public:
    /*! Construct a recursator from a matrix.
        \param matrix The matrix to which the recursator refers.
	\param bound  Explicit bound declaration; must not be smaller than the numbers of rows and the number of columns;
	              must also be a power of 2.

        Constructor takes the entire matrix as sub-matrix.
        This allows to have different type for the matrix and the sub-matrix.
    **/
    explicit matrix_recursator(Matrix const& matrix,
			       size_type bound= 0
			       ) 
	: my_sub_matrix(constructor_helper(matrix)), my_bound(outer_bound(matrix)),
	  my_first_row(0), my_first_col(0)         // splitter(*this)
    {
      if (bound == 0)
	my_bound= outer_bound(matrix);
      else {
	MTL_THROW_IF(!is_power_of_2(bound), range_error("Bound must be a power of 2"));
	MTL_THROW_IF(bound < matrix.num_rows() || bound < matrix.num_cols(), 
		     range_error("Bound must not be smaller than matrix dimensions"));
	my_bound= bound;
      }
    }


private:

    template <typename SubMatrix>
    sub_matrix_type get_value_dispatch(const SubMatrix& matrix, 
				   size_type br, size_type er, size_type bc, size_type ec) const
    {
	return sub_matrix(my_sub_matrix, br, er, bc, ec);
    }

    template <typename SubMatrix>
    sub_matrix_type get_value_dispatch(transposed_view<SubMatrix> view, 
				       size_type br, size_type er, size_type bc, size_type ec) const
    {
	typedef typename sub_matrix_t<SubMatrix>::sub_matrix_type   ref_sub_type;
	typedef boost::shared_ptr<ref_sub_type>                     pointer_type;
	typedef typename transposed_view<SubMatrix>::other          ref_type;

	pointer_type p(new ref_sub_type(sub_matrix(const_cast<ref_type&>(view.ref), bc, ec, br, er)));
	return sub_matrix_type(p); 	
    }


public:
    sub_matrix_type get_value() const
    {
	using std::min;
	size_type begin_row= my_sub_matrix.begin_row() + my_first_row,
	          end_row= min(begin_row + my_bound, my_sub_matrix.end_row()),
	          begin_col= my_sub_matrix.begin_col() + my_first_col,
	          end_col= min(begin_col + my_bound, my_sub_matrix.end_col());

#if 0
	std::cout << "get_value [" << begin_row << "-" << end_row << "]["
		  << begin_col << "-" << end_col << "]\n";
#endif
	return get_value_dispatch(my_sub_matrix, begin_row, end_row, begin_col, end_col);
    }

    /// Compute the sub-matrix corresponding to this recursator.
    sub_matrix_type operator*() const
    {
	return get_value();
    }

    // Returning quadrants for non-const recursator

    self north_west() const
    {
	self tmp(*this);
	tmp.my_bound >>= 1; // divide by 2
	return tmp;
    }

    self south_west() const
    {
	self tmp(*this);
	tmp.my_bound >>= 1; // divide by 2
	tmp.my_first_row += tmp.my_bound;
	return tmp;
    }

    self north_east() const
    {
	self tmp(*this);
	tmp.my_bound >>= 1; // divide by 2
	tmp.my_first_col += tmp.my_bound;
	return tmp;
    }

    self south_east() const
    {
	self tmp(*this);
	tmp.my_bound >>= 1; // divide by 20
	tmp.my_first_row += tmp.my_bound;
	tmp.my_first_col += tmp.my_bound;
	return tmp;
    }

    bool is_empty() const
    {
	return my_first_row >= my_sub_matrix.num_rows() || my_first_col >= my_sub_matrix.num_cols();
    }

    size_type num_rows() const
    {
	using std::min;
	size_type tmp= ::mtl::num_rows(my_sub_matrix);
	return my_first_row >= tmp ? 0 : min(my_bound, tmp - my_first_row);
    }
	
    size_type num_cols() const
    {
	using std::min;
	size_type tmp= ::mtl::num_cols(my_sub_matrix);
	return my_first_col >= tmp ? 0 : min(my_bound, tmp - my_first_col);
    }


    /// Return the bound of the recursator
    size_type bound() const
    {
	return my_bound;
    }

    /*! Set the bound of the recursator.
	\param bound  The new virtual bound; must be a power of 2.

        This function allows to declare a virtual bound smaller than the number of rows and/or columns.
	It must be used with uttermost care.
    **/
    void set_bound(size_type b)
    {
	my_bound= b;
    }

    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; /// Referred matrix (from which the sub-matrices are built)
    size_type           my_bound,      /// Virtual matrix size, i.e. upper bound for size of sub-matrix.
	                my_first_row,  /// Row of first entry in submatrix
                        my_first_col;  /// Row of first entry in submatrix 
};



// To use matrix_recursator with const matrices Reference must be 'Matrix const&'
template <typename Matrix, typename Splitter = max_dim_splitter<Matrix> >
struct matrix_recursator_s
{
    typedef matrix_recursator_s                                    self;
    typedef Matrix                                                matrix_type;
    typedef Splitter                                              splitter_type;
    typedef typename sub_matrix_t<Matrix>::sub_matrix_type        sub_matrix_type;
    typedef typename sub_matrix_t<Matrix>::const_sub_matrix_type  const_sub_matrix_type;
    typedef typename Matrix::size_type                            size_type;
    // typedef outer_bound_splitter<self>                            splitter_type;

private:
    
    // template <typename Matrix> why was it templated ???
    sub_matrix_type constructor_helper(Matrix const& matrix)
    {
	return sub_matrix(matrix, matrix.begin_row(), matrix.end_row(),
			  matrix.begin_col(), matrix.end_col());
    }

    // For views without own data, we need to generate a new sub_matrix as shared_ptr
    // template <typename Matrix>
    sub_matrix_type constructor_helper(transposed_view<Matrix> const& matrix)
    {
	typedef typename sub_matrix_t<Matrix>::sub_matrix_type   ref_sub_type;
	typedef boost::shared_ptr<ref_sub_type>                  pointer_type;

	// Submatrix of referred matrix, colums and rows interchanged
	// Create a submatrix, whos address will be kept by transposed_view
	pointer_type p(new ref_sub_type(sub_matrix(matrix.ref, matrix.begin_col(), matrix.end_col(), 
						   matrix.begin_row(), matrix.end_row())));
	return sub_matrix_type(p); 
    }

public:
    // Constructor takes the whole matrix as sub-matrix
    // This allows to have different type for the matrix and the sub-matrix
    // This also enables matrices to have references as sub-matrices
    explicit matrix_recursator_s(Matrix const& matrix, size_type bound= 0)