morton_dense.hpp

矩阵运算源码最新版本

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

    // Construction from difference of matrices
    template <typename E1, typename E2>
    morton_dense(const matrix::mat_mat_minus_expr<E1, E2>& src) : expr_base(*this)
    {
	change_dim(num_rows(src.first), num_cols(src.first));
	*this= src.first;
	*this-= src.second;
    }

    // Construction from product of matrices
    template <typename E1, typename E2>
    morton_dense(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)
    {
	MTL_THROW_IF(this->extern_memory && (num_rows != this->num_rows() || num_cols != this->num_cols()),
		     runtime_error("Can't change the size of matrices with external memory"));

	set_ranges(num_rows, num_cols);
	this->realloc(memory_need(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;
    }

    using assign_base::operator=;


    value_type operator() (key_type const& key) const
    {
	return this->data[key.dilated_row.dilated_value() + key.dilated_col.dilated_value()];
    }

    void operator()(key_type const& key, value_type const& value)
    {
	this->data[key.dilated_row.dilated_value() + key.dilated_col.dilated_value()]= value;
    }

    const_access_type operator() (size_type row, size_type col) const
    {
	return this->data[dilated_row_t(row).dilated_value() + dilated_col_t(col).dilated_value()];
    }

    value_type& operator() (size_type row, size_type col)
    {
	return this->data[dilated_row_t(row).dilated_value() + dilated_col_t(col).dilated_value()];
    }


  protected:
    void set_nnz()
    {
      this->my_nnz = this->num_rows() * this->num_cols();
    }
    
    size_type memory_need(size_type rows, size_type cols)
    {
        dilated_row_t n_rows(rows - 1);
        dilated_col_t n_cols(cols - 1);
        return (n_rows.dilated_value() + n_cols.dilated_value() + 1);
    }

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

    template <typename> friend struct sub_matrix_t;    

#if 0
    size_type my_used_memory;
    
  private:
  // add morton member variables here

    int rows_;          // number of rows of the matrix
    int cols_;          // number of columns of the matrix
    int quadOrder_;     // order of the level-0 quad
                        // quadOrder_ = pow(2, (int)log2(max(rows_, cols_))) + 1)
                        // or quadOrder = pow(2, (int)log2(max(rows_, cols_)))),
                        // depending on the value of rows_ and cols_.
    int storageSize_;   // size of allocated storage for the matrix
    int maxLevel_;      // maximum level of the quadtree for the matrix

    std::vector<int> upBoundVec_;    // upper boundary vector
    std::vector<int> lowBoundVec_;   // lower boundary vector
    std::vector<int> rowMaskVec_;    // row mask vector
    std::vector<int> colMaskVec_;    // col mask vector
    // T* data_;                   // a pointer to the matrix data array

    void setQuadOrder();        // set quadOrder_
    void setStorageSize();      // set storageSize_
    void setMaxLevel();         // set maxLevel_
    void setBoundVec();         // set boundary vectors
    void setMaskVec();          // set mask vectors
    void mkMortonSPDMatrix();   // make default Morton matrix
#endif
};


#if 0
// boundary checking of Ahnentafel index

// check if the index is the root
template <typename Elt, unsigned long BitMask, typename Parameters>
bool morton_dense<Elt, BitMask, Parameters>::isRoot(const AhnenIndex& index) const {
    return index.getIndex() == 3;
}

// check if the index is a leaf
template <typename Elt, unsigned long BitMask, typename Parameters>
bool morton_dense<Elt, BitMask, Parameters>::isLeaf(const AhnenIndex& index) const {
  // a possible better way: compare index with the boundary
  // vector directly, instead of calling isInBound()
    if(isInBound(index))
        return (index.getIndex() >= lowBoundVec_[maxLevel_]);
    else return 0;
}
#endif



// ================
// Free functions
// ================

template <typename Value, unsigned long Mask, typename Parameters>
typename morton_dense<Value, Mask, Parameters>::size_type
inline num_rows(const morton_dense<Value, Mask, Parameters>& matrix)
{
    return matrix.num_rows();
}

template <typename Value, unsigned long Mask, typename Parameters>
typename morton_dense<Value, Mask, Parameters>::size_type
inline num_cols(const morton_dense<Value, Mask, Parameters>& matrix)
{
    return matrix.num_cols();
}

template <typename Value, unsigned long Mask, typename Parameters>
typename morton_dense<Value, Mask, Parameters>::size_type
inline size(const morton_dense<Value, Mask, Parameters>& matrix)
{
    return matrix.num_cols() * matrix.num_rows();
}




// ================
// Range generators
// ================

namespace traits
{
    // VC 8.0 finds ambiguity with mtl::tag::morton_dense (I wonder why)
    using mtl::morton_dense;

    // ===========
    // For cursors
    // ===========

    template <class Elt, unsigned long BitMask, class Parameters>
    struct range_generator<glas::tag::all, morton_dense<Elt, BitMask, Parameters> >
    {
	typedef morton_dense<Elt, BitMask, Parameters>        Matrix;
	typedef complexity_classes::linear_cached        complexity;
	static int const                         level = 1;
	typedef morton_dense_el_cursor<BitMask>  type;
	type begin(Matrix const& matrix)
	{
	    return type(matrix.begin_row(), matrix.begin_col(), matrix.num_cols());
	}
	type end(Matrix const& matrix)
	{
	    return type(matrix.end_row(), matrix.begin_col(), matrix.num_cols());
	}
    };

    template <class Elt, unsigned long BitMask, class Parameters>
    struct range_generator<glas::tag::nz, morton_dense<Elt, BitMask, Parameters> >
	: range_generator<glas::tag::all, morton_dense<Elt, BitMask, Parameters> >
    {};

    template <class Elt, unsigned long BitMask, class Parameters>
    struct range_generator<glas::tag::row, morton_dense<Elt, BitMask, Parameters> >
	: detail::all_rows_range_generator<morton_dense<Elt, BitMask, Parameters>, complexity_classes::linear_cached> 
    {};

    // For a cursor pointing to some row give the range of elements in this row 
    template <class Elt, unsigned long BitMask, class Parameters>
    struct range_generator<glas::tag::nz, 
			   detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::row, 2> >
    {
	typedef morton_dense<Elt, BitMask, Parameters>                   matrix;
	typedef detail::sub_matrix_cursor<matrix, glas::tag::row, 2>  cursor;
	typedef complexity_classes::linear_cached                        complexity;
	static int const                                                 level = 1;
	typedef morton_dense_col_cursor<BitMask>                         type;
	
	type begin(cursor const& c)
	{
	    return type(c.key, c.ref.begin_col());
	}
	type end(cursor const& c)
	{
	    return type(c.key, c.ref.end_col());
	}
    };

    template <class Elt, unsigned long BitMask, class Parameters>
    struct range_generator<glas::tag::all, 
			   detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::row, 2> >
        : range_generator<glas::tag::nz, 
			  detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::row, 2> >
    {};

    template <class Elt, unsigned long BitMask, class Parameters>
    struct range_generator<glas::tag::col, morton_dense<Elt, BitMask, Parameters> >
	: detail::all_cols_range_generator<morton_dense<Elt, BitMask, Parameters>, complexity_classes::linear_cached> 
    {};

    // For a cursor pointing to some row give the range of elements in this row 
    template <class Elt, unsigned long BitMask, class Parameters>
    struct range_generator<glas::tag::nz, 
			   detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::col, 2> >
    {
	typedef morton_dense<Elt, BitMask, Parameters>                   matrix;
	typedef detail::sub_matrix_cursor<matrix, glas::tag::col, 2>  cursor;
	typedef complexity_classes::linear_cached                        complexity;
	static int const                                                 level = 1;
	typedef morton_dense_row_cursor<BitMask>                         type;
	
	type begin(cursor const& c)
	{
	    return type(c.ref.begin_row(), c.key);
	}
	type end(cursor const& c)
	{
	    return type(c.ref.end_row(), c.key);
	}
    };

    template <class Elt, unsigned long BitMask, class Parameters>
    struct range_generator<glas::tag::all, 
			   detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::col, 2> >
        : range_generator<glas::tag::nz, 
			  detail::sub_matrix_cursor<morton_dense<Elt, BitMask, Parameters>, glas::tag::col, 2> >
    {};


// =============
// For iterators
// =============

    namespace detail {

        template <typename OuterTag, typename Matrix, bool is_const>
        struct morton_dense_iterator_range_generator
        {
	    typedef Matrix                                                                matrix_type;
	    typedef typename matrix_type::size_type                                       size_type;
	    typedef typename matrix_type::value_type                                      value_type;
	    typedef typename matrix_type::parameters                                      parameters;
	    typedef detail::sub_matrix_cursor<matrix_type, OuterTag, 2>                   cursor;

	    typedef complexity_classes::linear_cached                                     complexity;
	    static int const                                                              level = 1;

	    typedef typename boost::mpl::if_<
		boost::is_same<OuterTag, glas::tag::row>
	      , typename boost::mpl::if_c<
    	            is_const 
		  , morton_dense_col_const_iterator<Matrix>
		  , morton_dense_col_iterator<Matrix>
		>::type
	      , typename boost::mpl::if_c<
    	            is_const 
		  , morton_dense_row_const_iterator<Matrix>
		  , morton_dense_row_iterator<Matrix>
		>::type
	    >::type type;  

        private:

	    typedef typename boost::mpl::if_c<is_const, const Matrix&, Matrix&>::type    mref_type; 

	    type begin_dispatch(cursor const& c, glas::tag::row)
	    {
		return type(const_cast<mref_type>(c.ref), c.key, c.ref.begin_col());
	    }
	    
	    type end_dispatch(cursor const& c, glas::tag::row)
	    {
		return type(const_cast<mref_type>(c.ref), c.key, c.ref.end_col());
	    }

	    type begin_dispatch(cursor const& c, glas::tag::col)
	    {
		return type(const_cast<mref_type>(c.ref), c.ref.begin_row(), c.key);
	    }

	    type end_dispatch(cursor const& c, glas::tag::col)
	    {
		return type(const_cast<mref_type>(c.ref), c.ref.end_row(), c.key);
	    }

        public:

	    type begin(cursor const& c)
	    {
		return begin_dispatch(c, OuterTag());
	    }

	    type end(cursor const& c)
	    {
		return end_dispatch(c, OuterTag());
	    }	
        };

    } // namespace detail

        
    template <typename Value, unsigned long BitMask, typename Parameters, typename OuterTag>
    struct range_generator<tag::iter::nz, 
			   detail::sub_matrix_cursor<morton_dense<Value, BitMask, Parameters>, OuterTag, 2> >
      : public detail::morton_dense_iterator_range_generator<OuterTag, morton_dense<Value, BitMask, Parameters>, false>
    {};

    template <typename Value, unsigned long BitMask, typename Parameters, typename OuterTag>
    struct range_generator<tag::iter::all, 
			   detail::sub_matrix_cursor<morton_dense<Value, BitMask, Parameters>, OuterTag, 2> >
      : public detail::morton_dense_iterator_range_generator<OuterTag, morton_dense<Value, BitMask, Parameters>, false>
    {};

    template <typename Value, unsigned long BitMask, typename Parameters, typename OuterTag>
    struct range_generator<tag::const_iter::nz, 
			   detail::sub_matrix_cursor<morton_dense<Value, BitMask, Parameters>, OuterTag, 2> >
      : public detail::morton_dense_iterator_range_generator<OuterTag, morton_dense<Value, BitMask, Parameters>, true>
    {};

    template <typename Value, unsigned long BitMask, typename Parameters, typename OuterTag>
    struct range_generator<tag::const_iter::all, 
			   detail::sub_matrix_cursor<morton_dense<Value, BitMask, Parameters>, OuterTag, 2> >
      : public detail::morton_dense_iterator_range_generator<OuterTag, morton_dense<Value, BitMask, Parameters>, true>
    {};


} // namespace traits


// ==========
// Sub matrix
// ==========

template <typename Value, unsigned long BitMask, typename Parameters>
struct sub_matrix_t<morton_dense<Value, BitMask, Parameters> >
{
    typedef morton_dense<Value, BitMask, Parameters>    matrix_type;
    typedef matrix_type                     sub_matrix_type;
    typedef matrix_type const               const_sub_matrix_type;
    typedef typename matrix_type::size_type size_type;
    
    sub_matrix_type operator()(matrix_type& matrix, size_type begin_r, size_type end_r, size_type begin_c, size_type end_c)
    {
	matrix.check_ranges(begin_r, end_r, begin_c, end_c);

	// Treat empty sub-matrices first (don't hold the memory contiguousness check (but don't need to))
	if (begin_r >= end_r || begin_c >= end_c) {
	    sub_matrix_type  tmp(matrix);
	    tmp.set_ranges(0, 0);
	    tmp.extern_memory= true;
	    return tmp;
	}

	// Check whether sub-matrix is contigous memory block
	// by comparing the address of the last and the first element in the entire and the sub-matrix
	MTL_DEBUG_THROW_IF(&matrix[end_r-1][end_c-1] - &matrix[begin_r][begin_c] 
			   != &matrix[end_r-begin_r-1][end_c-begin_c-1] - &matrix[0][0],
			   range_error("This sub-matrix cannot be used because it is split in memory"));
	// Check with David if this is a sufficient condition (it is a necessary at least)

	sub_matrix_type  tmp(matrix);

	typename matrix_type::dilated_row_t  dilated_row(begin_r);
	typename matrix_type::dilated_col_t  dilated_col(begin_c);

	// Set new start address within masked matrix
	tmp.data += dilated_row.dilated_value() + dilated_col.dilated_value();
	tmp.set_ranges(end_r - begin_r, end_c - begin_c);

	// sub matrix doesn't own the memory (and must not free at the end)
	tmp.extern_memory= true;

	return tmp;

    }

    const_sub_matrix_type
    operator()(matrix_type const& matrix, size_type begin_r, size_type end_r, size_type begin_c, size_type end_c)
    {
	// To minimize code duplication, we use the non-const version
	sub_matrix_type tmp((*this)(const_cast<matrix_type&>(matrix), begin_r, end_r, begin_c, end_c));
	return tmp;
    }	
};

} // namespace mtl

#endif // MTL_MORTON_DENSE_INCLUDE