dense2d.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_DENSE2D_INCLUDE
#define MTL_DENSE2D_INCLUDE

#include <algorithm>
#include <boost/type_traits.hpp>
#include <boost/mpl/if.hpp>

#include <boost/numeric/mtl/mtl_fwd.hpp>
#include <boost/numeric/mtl/utility/common_include.hpp>
#include <boost/numeric/mtl/detail/base_sub_matrix.hpp>
#include <boost/numeric/mtl/detail/contiguous_memory_block.hpp>
#include <boost/numeric/mtl/operation/set_to_zero.hpp>
#include <boost/numeric/mtl/utility/dense_el_cursor.hpp>
#include <boost/numeric/mtl/utility/strided_dense_el_cursor.hpp>
#include <boost/numeric/mtl/utility/strided_dense_el_iterator.hpp>
#include <boost/numeric/mtl/matrix/all_mat_expr.hpp>
#include <boost/numeric/mtl/matrix/operators.hpp>
#include <boost/numeric/mtl/operation/compute_factors.hpp>


namespace mtl {

using std::size_t;

// Forward declarations
template <typename Value, typename Parameters> class dense2D;
struct dense2D_indexer;


// Indexing for dense matrices
struct dense2D_indexer 
{
  private:
    // helpers for public functions
    size_t offset(size_t ldim, size_t r, size_t c, row_major) const 
    {
	return r * ldim + c; 
    }
    size_t offset(size_t ldim, size_t r, size_t c, col_major) const 
    {
	return c * ldim + r; 
    }
    
    size_t row(size_t offset, size_t ldim, row_major) const 
    {
	return offset / ldim; 
    }
    size_t row(size_t offset, size_t ldim, col_major) const 
    {
	return offset % ldim;
    }
    
    size_t col(size_t offset, size_t ldim, row_major) const 
    {
	return offset % ldim;
    }
    size_t col(size_t offset, size_t ldim, col_major) const 
    {
	return offset / ldim; 
    }

 public:
    template <typename Value, class Parameters>
    size_t operator() (const dense2D<Value, Parameters>& ma, size_t r, size_t c) const
    {
	typedef dense2D<Value, Parameters> matrix_type;
	// convert into c indices
	typename matrix_type::index_type my_index;
	size_t my_r= index::change_from(my_index, r);
	size_t my_c= index::change_from(my_index, c);
	return offset(ma.ldim, my_r, my_c, typename matrix_type::orientation());
    }

    template <typename Value, class Parameters>
    size_t row(const dense2D<Value, Parameters>& ma, 
	       typename dense2D<Value, Parameters>::key_type key) const
    {
	typedef dense2D<Value, Parameters> matrix_type;
	// row with c-index for my orientation
	size_t r= row(ma.offset(key), ma.ldim, typename matrix_type::orientation());
	return index::change_to(typename matrix_type::index_type(), r);
    }

    template <typename Value, class Parameters>
    size_t col(const dense2D<Value, Parameters>& ma, 
	       typename dense2D<Value, Parameters>::key_type key) const 
    {
	typedef dense2D<Value, Parameters> matrix_type;
	// column with c-index for my orientation
	size_t c= col(ma.offset(key), ma.ldim, typename matrix_type::orientation());
	return index::change_to(typename matrix_type::index_type(), c);
    }
    template <typename, typename> friend struct dense2D;
}; // dense2D_indexer


namespace detail 
{
    
    // Compute required memory
    // Enabling mechanism to make sure that computation is valid
    template <typename Parameters, bool Enable>
    struct dense2D_array_size {
	static std::size_t const value= 0;
    };

    template <typename Parameters>
    struct dense2D_array_size<Parameters, true>
    {
	typedef typename Parameters::dimensions   dimensions;
	BOOST_STATIC_ASSERT((dimensions::is_static));
	static std::size_t const value= dimensions::Num_Rows * dimensions::Num_Cols;
    };

} // namespace detail


// Forward declaration (for friend declaration)
namespace traits { namespace detail {
    template <typename, typename, bool> struct dense2D_iterator_range_generator;
}}

  
// Dense 2D matrix type
template <typename Value, typename Parameters = mtl::matrix::parameters<> >
class dense2D : public detail::base_sub_matrix<Value, Parameters>, 
		public detail::contiguous_memory_block< Value, Parameters::on_stack, 
							 detail::dense2D_array_size<Parameters, Parameters::on_stack>::value >,
                public detail::crtp_base_matrix< dense2D<Value, Parameters>, Value, std::size_t >,
		public matrix::mat_expr< dense2D<Value, Parameters> >
{
    typedef dense2D                                           self;
    typedef detail::base_sub_matrix<Value, Parameters>        super;
    typedef detail::contiguous_memory_block<Value, Parameters::on_stack, 
					     detail::dense2D_array_size<Parameters, Parameters::on_stack>::value>     super_memory;
    typedef matrix::mat_expr< dense2D<Value, Parameters> >    expr_base;
    typedef detail::crtp_base_matrix< self, Value, std::size_t > crtp_base;
    typedef detail::crtp_matrix_assign< self, Value, std::size_t > assign_base;
  public:
    typedef Parameters                        parameters;
    typedef typename Parameters::orientation  orientation;
    typedef typename Parameters::index        index_type;
    typedef typename Parameters::dimensions   dim_type;
    typedef Value                             value_type;
    // return type of operator() const
    typedef const value_type&                 const_access_type;

    typedef const value_type*                 const_pointer_type;
    typedef const_pointer_type                key_type;
    typedef std::size_t                       size_type;
    typedef dense_el_cursor<Value>            el_cursor_type;  
    typedef dense2D_indexer                   indexer_type;

    // Self-similar type unless dimension is fixed
    // Not supported for the moment
    typedef self                              sub_matrix_type;  

  protected:
    // Obviously, the next 3 functions must be called after setting dimensions
    void set_nnz()
    {
	this->my_nnz = this->num_rows() * this->num_cols();
    }

    void set_ldim(row_major)
    {
	ldim= this->num_cols();
    }

    void set_ldim(col_major)
    {
	ldim= this->num_rows();
    }

    void set_ldim()
    {
	set_ldim(orientation());
    }

    void init()
    {
      set_nnz(); set_ldim(); // set_to_zero(*this);
    }

  public:
    // if compile time matrix size allocate memory
    dense2D() : super(), super_memory(dim_type().num_rows() * dim_type().num_cols()), expr_base(*this) 
    { 
	init(); 
    }

    // only sets dimensions, only for run-time dimensions
    explicit dense2D(mtl::non_fixed::dimensions d) 
	: super(d), super_memory(d.num_rows() * d.num_cols()), expr_base(*this) 
    { 
	init(); 
    }

    dense2D(size_type num_rows, size_type num_cols) 
	: super(mtl::non_fixed::dimensions(num_rows, num_cols)), 
	  super_memory(num_rows * num_cols), expr_base(*this) 
    { 
	init(); 
    }

    // sets dimensions and pointer to external data
    explicit dense2D(mtl::non_fixed::dimensions d, value_type* a) 
      : super(d), super_memory(a, d.num_rows() * d.num_cols()), expr_base(*this) 
    { 
	init(); 
    }

    // sets dimensions and pointer to external data
    dense2D(size_type num_rows, size_type num_cols, value_type* a) 
      : super(mtl::non_fixed::dimensions(num_rows, num_cols)), super_memory(a, num_rows * num_cols), expr_base(*this) 
    { 
	init(); 
    }

    // same constructor for compile time matrix size
    // sets dimensions and pointer to external data
    explicit dense2D(value_type* a) 
	: super(), super_memory(a, dim_type().num_rows() * dim_type().num_cols()), expr_base(*this) 
    { 
	BOOST_STATIC_ASSERT((dim_type::is_static));
        init();
    }

    // Default copy constructor doesn't work because CRTP refers to copied matrix not to itself 
    dense2D(const self& m) 
	: super(mtl::non_fixed::dimensions(m.num_rows(), m.num_cols())), 
	  super_memory(&(const_cast<self&>(m)[0][0]), m.num_rows() * m.num_cols()), expr_base(*this)
    {
	//init();
	this->my_nnz= m.my_nnz; ldim= m.ldim;
	// std::cout << "In copy constructor:\n"; print_matrix(*this);
    }

#ifndef _MSC_VER // Constructors need rigorous reimplementation, cf. #142-#144
    // Construction from sum of matrices
    template <typename E1, typename E2>
    dense2D(const matrix::mat_mat_plus_expr<E1, E2>& src)
	: super(mtl::non_fixed::dimensions(num_rows(src.first), num_cols(src.first))),
	  super_memory(num_rows(src.first) * num_cols(src.first)), expr_base(*this)		
    {
	init();
	matrix_copy(src.first, *this);
	*this+= src.second;
    }

    // Construction from difference of matrices
    template <typename E1, typename E2>
    dense2D(const matrix::mat_mat_minus_expr<E1, E2>& src)
	: super(mtl::non_fixed::dimensions(num_rows(src.first), num_cols(src.first))),
	  super_memory(num_rows(src.first) * num_cols(src.first)), expr_base(*this)		
    {
	init();
	matrix_copy(src.first, *this);
	*this-= src.second;
    }


    // Construction from product of matrices
    template <typename E1, typename E2>
    dense2D(const matrix::mat_mat_times_expr<E1, E2>& src)
	: expr_base(*this)		
    {
	operation::compute_factors<self, matrix::mat_mat_times_expr<E1, E2> > factors(src);
	change_dim(num_rows(factors.first), num_cols(factors.second));
	mult(factors.first, factors.second, *this);
    }
#endif


    void change_dim(size_type num_rows, size_type num_cols)
    {
	super::change_dim(mtl::non_fixed::dimensions(num_rows, num_cols));
	set_nnz(); set_ldim();
	this->realloc(num_rows * num_cols);
    }


    // Alleged ambiguity in MSVC 8.0, I need to turn off the warning 
	// Removing the operator ends in run-time error
   self& operator=(const self& src)
    {
	// no self-copy
	if (this == &src) return *this;

	change_dim(src.num_rows(), src.num_cols());
	std::copy(src.elements(), src.elements()+src.used_memory(), this->elements());
	return *this;
    }

    // import operators from CRTP base class
    using assign_base::operator=;


    bool check_indices(size_t r, size_t c) const
    {
	return r >= this->begin_row() && r < this->end_row() && c >= this->begin_col() && c < this->end_col();
    }

    
    const_access_type operator() (size_t r, size_t c) const 
    {
	// assert(check_indices(r, c));  // causes trouble for iterator/cursor creation
	size_t offset= indexer(*this, r, c);
        return this->data[offset];
    }

    value_type& operator() (size_t r, size_t c)
    {
	// assert(check_indices(r, c));  // causes trouble for iterator/cursor creation
	return this->data[indexer(*this, r, c)]; 
    }    

    // offset regarding c-style indices
    size_t c_offset(size_t r, size_t c) const
    {
	return indexer.offset(ldim, r, c, orientation());
    }

    size_type get_ldim() const
    {
	return ldim;
    }

    friend void swap(self& matrix1, self& matrix2)
    {
	static_cast<super_memory&>(matrix1).swap(matrix2);
	static_cast<super&>(matrix1).swap(matrix2);
	std::swap(matrix1.ldim, matrix2.ldim);
    }


  protected:
    
    // Set ranges from begin_r to end_r and begin_c to end_c
    void set_ranges(size_type begin_r, size_type end_r, size_type begin_c, size_type end_c)
    {
	super::set_ranges(begin_r, end_r, begin_c, end_c);
	set_nnz();
    }
	
    // Set ranges to a num_row x num_col matrix, keeps indexing
    void set_ranges(size_type num_rows, size_type num_cols)
    {
	set_ranges(this->begin_row(), this->begin_row() + num_rows, 
		   this->begin_col(), this->begin_col() + num_cols);
    }
    

  public:

    indexer_type  indexer;
  protected:
    // Leading dimension is minor dimension in original matrix 
    // Opposed to other dims doesn't change in sub-matrices
    size_type     ldim; 

    friend class dense2D_indexer;
    template <typename> friend struct sub_matrix_t;
    template <typename, typename> friend struct traits::range_generator;
    template <typename, typename, bool> friend struct traits::detail::dense2D_iterator_range_generator;