compressed2d.hpp

矩阵运算源码最新版本

    typedef size_t                                   size_type;
    typedef compressed2D_indexer                     indexer_type;

    // Only allocation of new data, doesn't copy if already existent
    void allocate(size_t new_nnz)
    {
	if (new_nnz) {
	    this->my_nnz = new_nnz;
	    data.resize(this->my_nnz);
	    indices.resize(this->my_nnz, 0);
	}
    }

    // removes all values; e.g. for set_to_zero
    void make_empty()
    {
	this->my_nnz = 0;
	data.resize(0);
	indices.resize(0);
	std::fill(starts.begin(), starts.end(), 0);
    }

    void change_dim(size_type num_rows, size_type num_cols)
    {
	super::change_dim(mtl::non_fixed::dimensions(num_rows, num_cols));
	starts.resize(this->dim1()+1);
	make_empty();
    }

    // if compile time matrix size, we can set the start vector
    explicit compressed2D () 
	: super(), expr_base(*this), inserting(false)
    {
	if (super::dim_type::is_static) starts.resize(super::dim1() + 1);
    }

    // setting dimension and allocate starting vector
    explicit compressed2D (mtl::non_fixed::dimensions d, size_t nnz = 0) 
      : super(d), expr_base(*this), inserting(false)
    {
	starts.resize(super::dim1() + 1, 0);
	allocate(nnz);
    }

    // setting dimension and allocate starting vector
    compressed2D (size_type num_rows, size_type num_cols, size_t nnz = 0) 
      : super(non_fixed::dimensions(num_rows, num_cols)), expr_base(*this), inserting(false)
    {
	starts.resize(super::dim1() + 1, 0);
	allocate(nnz);
    }

    explicit compressed2D(const self& src)
      : super(non_fixed::dimensions(::mtl::num_rows(src), ::mtl::num_cols(src))), expr_base(*this), inserting(false)
    {
	starts.resize(super::dim1() + 1, 0);
	matrix_copy(src, *this);
    }

    template <typename SrcValue, typename SrcParameters>
    explicit compressed2D(const compressed2D<SrcValue, SrcParameters>& src)
	: super(non_fixed::dimensions(::mtl::num_rows(src), ::mtl::num_cols(src))), expr_base(*this), inserting(false)
    {
	starts.resize(super::dim1() + 1, 0);
	matrix_copy(src, *this);
    }


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

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

    // Construction from product of matrices
    template <typename E1, typename E2>
    explicit compressed2D(const matrix::mat_mat_times_expr<E1, E2>& src) 
	: expr_base(*this), inserting(false)		
    {
	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


    // 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;

	matrix_copy(src, *this);
	return *this;
    }

    using assign_base::operator=;


    // Copies range of values and their coordinates into compressed matrix
    // For brute force initialization, should be used with uttermost care
    // Won't be suitable for distributed matrices, take care of this to this later
    template <typename ValueIterator, typename StartIterator, typename IndexIterator>    
    void raw_copy(ValueIterator first_value, ValueIterator last_value, 
		  StartIterator first_start, IndexIterator first_index)
    {
	using std::copy;

	// check if starts has right size
	allocate(last_value - first_value); // ???? 
	// check if nnz and indices has right size

	copy(first_value, last_value, data.begin());
	copy(first_start, first_start + this->dim1() + 1, starts.begin());
	copy(first_index, first_index + this->nnz(), indices.begin());
    }

    // Consistency check urgently needed !!!

    const_access_type operator() (size_type row, size_type col) const
    {
        MTL_DEBUG_THROW_IF(inserting, logic_error("Reading data during insertion has undefined behavior"));
	utilities::maybe<size_type> pos = indexer(*this, row, col);
	return pos ? data[pos] : value_type(0);
    }

    value_type value_from_offset(size_type offset) const
    {
	MTL_DEBUG_THROW_IF(offset >= this->my_nnz, index_out_of_range("Offset larger than matrix"));
	return data[offset];
    }

    value_type& value_from_offset(size_type offset)
    {
	MTL_DEBUG_THROW_IF(offset >= this->my_nnz, index_out_of_range("Offset larger than matrix"));
	return data[offset];
    }

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

	swap(matrix1.data, matrix2.data);
	swap(matrix1.starts, matrix2.starts);
	swap(matrix1.indices, matrix2.indices);
	swap(matrix1.inserting, matrix2.inserting);
    }

    friend struct compressed2D_indexer;
    template <typename, typename, typename> friend struct compressed2D_inserter;
    template <typename, typename> friend struct compressed_el_cursor;
    template <typename, typename> friend struct compressed_minor_cursor;

    indexer_type            indexer;
    std::vector<value_type> data; 
  protected:
    std::vector<size_t>     starts;
    std::vector<size_t>     indices;
    bool                    inserting;
};



// ========
// Inserter
// ========

// Additional data structure to insert into compressed 2D matrix type
template <typename Elt, typename Parameters, typename Updater = mtl::operations::update_store<Elt> >
struct compressed2D_inserter
{
    typedef compressed2D_inserter             self;
    typedef compressed2D<Elt, Parameters>     matrix_type;
    typedef typename matrix_type::size_type   size_type;
    typedef typename matrix_type::value_type  value_type;
    typedef std::pair<size_type, size_type>   size_pair;
    typedef std::map<size_pair, value_type>   map_type;
    typedef operations::update_proxy<self, size_type>   proxy_type;

  private: 
    // stretch matrix rows or columns to slot size (or leave it if equal or greater)
    void stretch();

  public:
    explicit compressed2D_inserter(matrix_type& matrix, size_type slot_size = 5)
	: matrix(matrix), elements(matrix.data), starts(matrix.starts), indices(matrix.indices), 
	  slot_size(slot_size), slot_ends(matrix.dim1()) 
    {
	MTL_THROW_IF(matrix.inserting, runtime_error("Two inserters on same matrix"));
	matrix.inserting = true;
	stretch();
    }

    ~compressed2D_inserter()
    {
	final_place();
	insert_spare();
	matrix.inserting = false;
    }
	
    proxy_type operator() (size_type row, size_type col)
    {
	return proxy_type(*this, row, col);
    }

    void update(size_type row, size_type col, value_type val);

    template <typename Matrix, typename Rows, typename Cols>
    self& operator<< (const matrix::element_matrix_t<Matrix, Rows, Cols>& elements)
    {
	for (unsigned ri= 0; ri < elements.rows.size(); ri++)
	    for (unsigned ci= 0; ci < elements.cols.size(); ci++)
		update (elements.rows[ri], elements.cols[ci], elements.matrix(ri, ci));
	return *this;
    }

    template <typename Matrix, typename Rows, typename Cols>
    self& operator<< (const matrix::element_array_t<Matrix, Rows, Cols>& elements)
    {
	for (unsigned ri= 0; ri < elements.rows.size(); ri++)
	    for (unsigned ci= 0; ci < elements.cols.size(); ci++)
		update (elements.rows[ri], elements.cols[ci], elements.array[ri][ci]);
	return *this;
    }

  private:
	utilities::maybe<typename self::size_type> matrix_offset(size_pair);
    void final_place();
    void insert_spare();

  protected:
    compressed2D<Elt, Parameters>&      matrix;
    std::vector<value_type>&            elements;
    std::vector<size_type>&             starts;
    std::vector<size_type>&             indices;
    size_type                           slot_size;
    std::vector<size_type>              slot_ends;
    map_type                            spare;
};

template <typename Elt, typename Parameters, typename Updater>
void compressed2D_inserter<Elt, Parameters, Updater>::stretch()
{
    using std::copy;
    using std::copy_backward;
    using std::swap;

    std::vector<size_type>  new_starts(matrix.dim1() + 1);
    new_starts[0] = 0;
    for (size_type i = 0; i < matrix.dim1(); i++) {
	size_type entries = starts[i+1] - starts[i];
	slot_ends[i] = new_starts[i] + entries; 
	new_starts[i+1] = new_starts[i] + std::max(entries, slot_size);
    }

    size_type new_total = new_starts[matrix.dim1()];
    elements.resize(new_total);
    indices.resize(new_total);
	if (elements.empty()) return;
   
    // copy normally if not overlapping and backward if overlapping
    // i goes down to 1 (not to 0) because i >= 0 never stops for unsigned ;-)
	// &v[i] is replaced by &v[0]+i to enable past-end addresses for STL copy
    for (size_type i = matrix.dim1(); i > 0; i--)
	if (starts[i] <= new_starts[i-1]) {
	    copy(&elements[0] + starts[i-1], &elements[0] + starts[i], &elements[0] + new_starts[i-1]);
	    copy(&indices[0] + starts[i-1], &indices[0] + starts[i], &indices[0] + new_starts[i-1]);
	} else {
		copy_backward(&elements[0] + starts[i-1], &elements[0] + starts[i], &elements[0] + slot_ends[i-1]);
		copy_backward(&indices[0] + starts[i-1], &indices[0] + starts[i], &indices[0] + slot_ends[i-1]);
	}
    swap(starts, new_starts);		    
}

template <typename Elt, typename Parameters, typename Updater>
inline utilities::maybe<typename compressed2D_inserter<Elt, Parameters, Updater>::size_type> 
compressed2D_inserter<Elt, Parameters, Updater>::matrix_offset(size_pair mm)
{
    size_type major, minor;
    boost::tie(major, minor) = mm;