sparse.hpp

Go to the documentation of this file.

/*!
 * \brief       Classes used for vector proxies
 * 
 * \author      O. Krause
 * \date        2016
 *
 *
 * \par Copyright 1995-2015 Shark Development Team
 * 
 * <BR><HR>
 * This file is part of Shark.
 * <http://image.diku.dk/shark/>
 * 
 * Shark is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published 
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * Shark is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with Shark.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 #ifndef REMORA_SPARSE_HPP
#define REMORA_SPARSE_HPP
 
 
#include "detail/storage.hpp"
#include "detail/traits.hpp"
#include "cpu/sparse.hpp"
#include "cpu/sparse_matrix.hpp"
#include "expression_types.hpp"
#include "assignment.hpp"
 
 
namespace remora{
 
template<class T, class I = std::size_t> class compressed_vector;
template<class T, class I = std::size_t> class compressed_vector_adaptor;
template<class T, class I = std::size_t, class Orientation = row_major> class compressed_matrix;
 
 
/** \brief Compressed array based sparse vector
 *
 * a sparse vector of values of type T of variable size. The non zero values are stored as
 * two seperate arrays: an index array and a value array. The index array is always sorted
 * and there is at most one entry for each index. Inserting an element can be time consuming.
 * If the vector has a very high dimension with a few non-zero values, then this vector is
 * very time and memory efficient.
 *
 * For a \f$n\f$-dimensional compressed vector and \f$0 \leq i < n\f$ the non-zero elements
 * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for
 * elements \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of these containers holds \f$i_1 < i_2\f$.
 *
 * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
 * \tparam I the indices stored in the vector
 */
template<class T, class I>
class compressed_vector
:public vector_container<compressed_vector<T, I>, cpu_tag >
,public detail::BaseSparseVector<detail::VectorStorage<T,I > >{
public:
    typedef T value_type;
    typedef I size_type;
    typedef T const& const_reference;
    typedef T& reference;
    
    typedef detail::compressed_vector_reference<compressed_vector const> const_closure_type;
    typedef detail::compressed_vector_reference<compressed_vector> closure_type;
    typedef sparse_vector_storage<T,I> storage_type;
    typedef sparse_vector_storage<T const,I const> const_storage_type;
    typedef elementwise<sparse_tag> evaluation_category;
 
    // Construction and destruction
    compressed_vector()
    : detail::BaseSparseVector<detail::VectorStorage<T,I> >(
        detail::VectorStorage<T, I>(),0,0
    ){}
    explicit compressed_vector(size_type size, size_type non_zeros = 0)
    : detail::BaseSparseVector<detail::VectorStorage<T,I> >(
        detail::VectorStorage<T, I>(),size,non_zeros
    ){}
    template<class E>
    compressed_vector(vector_expression<E, cpu_tag> const& e, size_type non_zeros = 0)
    : detail::BaseSparseVector<detail::VectorStorage<T,I> >(
        detail::VectorStorage<T, I>(),e().size(),non_zeros
    ){
        assign(*this,e);
    }
    
    void resize(size_type size, bool keep = false){
        this->do_resize(size, keep);
    }
    
    ///\brief Returns the underlying storage structure for low level access
    storage_type raw_storage(){
        return this->m_manager.m_storage;
    }
    
    ///\brief Returns the underlying storage structure for low level access
    const_storage_type raw_storage() const{
        return this->m_manager.m_storage;
    }
    
    typename device_traits<cpu_tag>::queue_type& queue() const{
        return device_traits<cpu_tag>::default_queue();
    }
 
    friend void swap(compressed_vector& v1, compressed_vector& v2){
        std::swap(v1.m_size, v2.m_size);
        v1.m_indices.swap(v2.m_indices);
        v1.m_values.swap(v2.m_values);
        v1.m_storage.swap(v2.m_storage);
    }
 
    // Assignment
    compressed_vector& operator = (compressed_vector const& v) = default;
    template<class C>          // Container assignment without temporary
    compressed_vector& operator = (vector_container<C, cpu_tag> const& v) {
        this->resize(v().size(), false);
        assign(*this, v);
        return *this;
    }
    template<class AE>
    compressed_vector& operator = (vector_expression<AE, cpu_tag> const& ae) {
        compressed_vector temporary(ae, this->nnz_capacity());
        swap(temporary);
        return *this;
    }
 
    // Serialization
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /* file_version */) {
        ar & static_cast< detail::BaseSparseVector<detail::VectorStorage<T,I > >& >(*this);
    }
};
 
 
/** \brief Wraps external memory compatible to the format of a compressed vector
 *
 * For a \f$n\f$-dimensional compressed vector and \f$0 \leq i < n\f$ the non-zero elements
 * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for
 * elements \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of these containers holds \f$i_1 < i_2\f$.
 *
 * There are 4 values needed: the address of the arrays of indices and values, the number of nonzero elements
 * and the size of the arrays (which can be larger than the number of noznero elements to allow for insertion).
 *
 * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
 * \tparam I the indices stored in the vector
 */
template<class T, class I>
class compressed_vector_adaptor
: public vector_expression<compressed_vector_adaptor<T, I>, cpu_tag>
,public detail::BaseSparseVector<detail::VectorStorageReference<T,I > >{
public:
    typedef typename std::remove_const<T>::type value_type;
    typedef typename std::remove_const<I>::type size_type;
    typedef T const& const_reference;
    typedef T& reference;
    
    typedef detail::compressed_vector_reference<compressed_vector_adaptor const> const_closure_type;
    typedef detail::compressed_vector_reference<compressed_vector_adaptor> closure_type;
    typedef sparse_vector_storage<T,I> storage_type;
    typedef sparse_vector_storage<T const,I const> const_storage_type;
    typedef elementwise<sparse_tag> evaluation_category;
 
    // Construction and destruction
 
    /// \brief Constructs the adaptor from external storage
    explicit compressed_vector_adaptor(size_type size, storage_type storage)
    : detail::BaseSparseVector<detail::VectorStorageReference<T,I> >(
        detail::VectorStorageReference<T,I>(storage), size, storage.nnz
    ){}
    
    /// \brief Covnerts an expression into an adaptor
    template<class E>
    compressed_vector_adaptor(vector_expression<E, cpu_tag> const& e)
    : detail::BaseSparseVector<detail::VectorStorageReference<T,I> >(
        detail::VectorStorageReference<T,I>(e().raw_storage()), e().size(), e().raw_storage().nnz
    ){}
    
    ///\brief Returns the underlying storage structure for low level access
    storage_type raw_storage(){
        return this->m_manager.m_storage;
    }
    
    ///\brief Returns the underlying storage structure for low level access
    const_storage_type raw_storage() const{
        return this->m_manager.m_storage;
    }
    
    typename device_traits<cpu_tag>::queue_type& queue() const{
        return device_traits<cpu_tag>::default_queue();
    }
 
    // Assignment
    compressed_vector_adaptor& operator = (compressed_vector_adaptor const& v) {
        kernels::assign(*this, typename vector_temporary<compressed_vector_adaptor>::type(v));
        return *this;
    }
    template<class AE>
    compressed_vector_adaptor& operator = (vector_expression<AE, cpu_tag> const& ae) {
        kernels::assign(*this, typename vector_temporary<AE>::type(ae));
        return *this;
    }
};
 
template<class T, class I, class Orientation>
class compressed_matrix:public matrix_container<compressed_matrix<T, I, Orientation>, cpu_tag >{
public:
    typedef I size_type;
    typedef T value_type;
    typedef T const& const_reference;
    typedef T& reference;
    
    typedef detail::compressed_matrix_proxy<compressed_matrix<T, I, Orientation> const, Orientation> const_closure_type;
    typedef detail::compressed_matrix_proxy<compressed_matrix<T, I, Orientation>, Orientation> closure_type;
    typedef sparse_matrix_storage<T, I> storage_type;
    typedef sparse_matrix_storage<T const, I const> const_storage_type;
    typedef elementwise<sparse_tag> evaluation_category;
    typedef Orientation orientation;
 
    compressed_matrix():m_impl(detail::MatrixStorage<T,I>(0,0)){}
    
    compressed_matrix(size_type rows, size_type cols, size_type non_zeros = 0)
    :m_impl(detail::MatrixStorage<T,I>(orientation::index_M(rows,cols),orientation::index_m(rows,cols)),non_zeros){}
    
    template<class E>
    compressed_matrix(matrix_expression<E, cpu_tag> const& m, size_type non_zeros = 0)
    :m_impl(
        detail::MatrixStorage<T,I>(orientation::index_M(m().size1(),m().size2()),orientation::index_m(m().size1(),m().size2())),
        non_zeros
    ){
        assign(*this,m);
    }
    
    template<class E>
    compressed_matrix operator=(matrix_container<E, cpu_tag> const& m){
        resize(m.size1(),m.size2());
        return assign(*this,m);
    }
    template<class E>
    compressed_matrix& operator=(matrix_expression<E, cpu_tag> const& m){
        compressed_matrix temporary(m);
        m_impl = std::move(temporary.m_impl);
        return *this;
    }
    
    ///\brief Number of rows of the matrix
    size_type size1() const {
        return orientation::index_M(m_impl.major_size(), m_impl.minor_size());
    }
    
    ///\brief Number of columns of the matrix
    size_type size2() const {
        return orientation::index_m(m_impl.major_size(), m_impl.minor_size());
    }
 
    /// \brief Number of nonzeros this matrix can maximally store before requiring new memory
    std::size_t nnz_capacity() const{
        return m_impl.nnz_capacity();
    }
    /// \brief Number of reserved elements in the matrix (> number of nonzeros stored)
    std::size_t nnz_reserved() const {
        return m_impl.nnz_reserved();
    }
    /// \brief Number of nonzeros the major index (a row or column depending on orientation) can maximally store before a resize
    std::size_t major_capacity(size_type i)const{
        return m_impl.major_capacity(i);
    }
    /// \brief Number of nonzeros the major index (a row or column depending on orientation) currently stores
    std::size_t major_nnz(size_type i) const {
        return m_impl.major_nnz(i);
    }
 
    /// \brief Set the total number of nonzeros stored by the matrix
    void set_nnz(std::size_t non_zeros) {
        m_impl.set_nnz(non_zeros);
    }
    /// \brief Set the number of nonzeros stored in the major index (a row or column depending on orientation)
    void set_major_nnz(size_type i,std::size_t non_zeros) {
        m_impl.set_major_nnz(i,non_zeros);
    }
    
    const_storage_type raw_storage()const{
        return m_impl.raw_storage();
    }
    storage_type raw_storage(){
        return m_impl.raw_storage();
    }
    
    typename device_traits<cpu_tag>::queue_type& queue() const{
        return device_traits<cpu_tag>::default_queue();
    }
    
    void reserve(std::size_t non_zeros) {
        m_impl.reserve(non_zeros);
    }
 
    void major_reserve(size_type i, std::size_t non_zeros, bool exact_size = false) {
        m_impl.major_reserve(i, non_zeros, exact_size);
    }
 
    void resize(size_type rows, size_type columns){
        m_impl.resize(orientation::index_M(rows,columns),orientation::index_m(rows,columns));
    }
    
    typedef typename detail::compressed_matrix_impl<detail::MatrixStorage<T,I> >::const_major_iterator const_major_iterator;
    typedef typename detail::compressed_matrix_impl<detail::MatrixStorage<T,I> >::major_iterator major_iterator;
 
    const_major_iterator major_begin(size_type i) const {
        return m_impl.cmajor_begin(i);
    }
 
    const_major_iterator major_end(size_type i) const{
        return m_impl.cmajor_end(i);
    }
 
    major_iterator major_begin(size_type i) {
        return m_impl.major_begin(i);
    }
 
    major_iterator major_end(size_type i) {
        return m_impl.major_end(i);
    }
    
    major_iterator set_element(major_iterator pos, size_type index, value_type value){
        return m_impl.set_element(pos, index, value);
    }
 
    major_iterator clear_range(major_iterator start, major_iterator end) {
        return m_impl.clear_range(start,end);
    }
    
    void clear() {
        for(std::size_t i = 0; i != m_impl.major_size(); ++i){
            clear_range(major_begin(i),major_end(i));
        }
    }
    // Serialization
    template<class Archive>
    void serialize(Archive &ar, const unsigned int /* file_version */) {
        ar & m_impl;
    }
 
private:
    detail::compressed_matrix_impl<detail::MatrixStorage<T,I> > m_impl;
};
 
 
 
//~ ///\brief Wraps externally provided storage into a sparse matrix interface
//~ ///
//~ /// Note that for this class, storage is limited and if an insertion operation takes more space than available, it will throw an exception
//~ template<class T, class I, class Orientation>
//~ class compressed_matrix_adaptor:public matrix_container<compressed_matrix_adaptor<T, I, Orientation>, cpu_tag >{
//~ public:
    //~ typedef I size_type;
    //~ typedef T value_type;
    //~ typedef T const& const_reference;
    //~ typedef T& reference;
    
    //~ typedef detail::compressed_matrix_proxy<compressed_matrix<T, I, Orientation> const, Orientation> const_closure_type;
    //~ typedef detail::compressed_matrix_proxy<compressed_matrix<T, I, Orientation>, Orientation> closure_type;
    //~ typedef sparse_matrix_storage<T, I> storage_type;
    //~ typedef sparse_matrix_storage<T const, I const> const_storage_type;
    //~ typedef elementwise<sparse_tag> evaluation_category;
    //~ typedef Orientation orientation;
    
    //~ compressed_matrix_adaptor(size_type rows, size_type cols, storage_type storage)
    //~ :m_impl(detail::MatrixStorageAdaptor<T,I>(orientation::index_M(rows,cols),orientation::index_m(rows,cols)),storage){}
    
    //~ template<class E>
    //~ compressed_matrix operator=(matrix_container<E, cpu_tag> const& m){
        //~ return assign(*this,m);
    //~ }
    //~ template<class E>
    //~ compressed_matrix& operator=(matrix_expression<E, cpu_tag> const& m){
        //~ compressed_matrix temporary(m);
        //~ swap(*this,temporary);
        //~ return *this;
    //~ }
    
    //~ ///\brief Number of rows of the matrix
    //~ size_type size1() const {
        //~ return orientation::index_M(m_impl.major_size(), m_impl.minor_size());
    //~ }
    
    //~ ///\brief Number of columns of the matrix
    //~ size_type size2() const {
        //~ return orientation::index_m(m_impl.major_size(), m_impl.minor_size());
    //~ }
 
    //~ /// \brief Number of nonzeros this matrix can maximally store before memory is exhausted
    //~ std::size_t nnz_capacity() const{
        //~ return m_impl.nnz_capacity();
    //~ }
    //~ /// \brief Number of reserved elements in the matrix (> number of nonzeros stored, < nnz_capacity)
    //~ std::size_t nnz_reserved() const {
        //~ return m_impl.nnz_reserved();
    //~ }
    //~ /// \brief Number of nonzeros the major index (a row or column depending on orientation) can maximally store before a resize
    //~ std::size_t major_capacity(size_type i)const{
        //~ return m_impl.major_capacity(i);
    //~ }
    //~ /// \brief Number of nonzeros the major index (a row or column depending on orientation) currently stores
    //~ std::size_t major_nnz(size_type i) const {
        //~ return m_impl.major_nnz(i);
    //~ }
 
    //~ /// \brief Set the total number of nonzeros stored by the matrix
    //~ void set_nnz(std::size_t non_zeros) {
        //~ m_impl.set_nnz(non_zeros);
    //~ }
    //~ /// \brief Set the number of nonzeros stored in the major index (a row or column depending on orientation)
    //~ void set_major_nnz(size_type i,std::size_t non_zeros) {
        //~ m_impl.set_major_nnz(i,non_zeros);
    //~ }
    
    //~ const_storage_type raw_storage()const{
        //~ return m_impl.raw_storage();
    //~ }
    //~ storage_type raw_storage(){
        //~ return m_impl.raw_storage();
    //~ }
    
    //~ typename device_traits<cpu_tag>::queue_type& queue() const{
        //~ return device_traits<cpu_tag>::default_queue();
    //~ }
    
    //~ void reserve(std::size_t non_zeros) {
        //~ m_impl.reserve(non_zeros);
    //~ }
 
    //~ void major_reserve(size_type i, std::size_t non_zeros) {
        //~ m_impl.major_reserve(i, non_zeros, true);
    //~ }
    
    //~ typedef typename detail::compressed_matrix_impl<detail::MatrixStorage<T,I> >::const_major_iterator const_major_iterator;
    //~ typedef typename detail::compressed_matrix_impl<detail::MatrixStorage<T,I> >::major_iterator major_iterator;
 
    //~ const_major_iterator major_begin(size_type i) const {
        //~ return m_impl.cmajor_begin(i);
    //~ }
 
    //~ const_major_iterator major_end(size_type i) const{
        //~ return m_impl.cmajor_end(i);
    //~ }
 
    //~ major_iterator major_begin(size_type i) {
        //~ return m_impl.major_begin(i);
    //~ }
 
    //~ major_iterator major_end(size_type i) {
        //~ return m_impl.major_end(i);
    //~ }
    
    //~ major_iterator set_element(major_iterator pos, size_type index, value_type value){
        //~ return m_impl.set_element(pos, index, value);
    //~ }
 
    //~ major_iterator clear_range(major_iterator start, major_iterator end) {
        //~ return m_impl.clear_range(start,end);
    //~ }
    
    //~ void clear() {
        //~ for(std::size_t i = 0; i != m_impl.major_size(); ++i){
            //~ clear_range(major_begin(i),major_end(i));
        //~ }
    //~ }
//~ private:
    //~ detail::compressed_matrix_impl<detail::MatrixStorage<T,I> > m_impl;
//~ };
 
namespace detail{
////////////////////////MATRIX ROW//////////////////////
template<class M>
struct matrix_row_optimizer<detail::compressed_matrix_proxy<M, row_major> >{
    typedef compressed_matrix_row<detail::compressed_matrix_proxy<M, row_major> > type;
    
    static type create(detail::compressed_matrix_proxy<M, row_major> const& m, std::size_t i){
        //create vector reference
        return type(m,i);
    }
};
 
 
////////////////////////MATRIX TRANSPOSE//////////////////////
template<class M, class Orientation>
struct matrix_transpose_optimizer<detail::compressed_matrix_proxy<M, Orientation> >{
    typedef detail::compressed_matrix_proxy<M, typename Orientation::transposed_orientation> type;
    
    static type create(detail::compressed_matrix_proxy<M, Orientation> const& m){
        return type(m.matrix(), m.start(), m.end());
    }
};
 
////////////////////////MATRIX ROWS//////////////////////
template<class M>
struct matrix_rows_optimizer<detail::compressed_matrix_proxy<M, row_major> >{
    typedef detail::compressed_matrix_proxy<M, row_major> type;
    
    static type create(detail::compressed_matrix_proxy<M, row_major> const& m, 
        std::size_t start, std::size_t end
    ){
        return type(m.matrix(), start + m.start(), end+m.start());
    }
};
}
 
template<class T, class O>
struct matrix_temporary_type<T,O,sparse_tag, cpu_tag> {
    typedef compressed_matrix<T,std::size_t, O> type;
};
 
template<class T>
struct vector_temporary_type<T,sparse_tag, cpu_tag>{
    typedef compressed_vector<T> type;
};
 
}
#endif