matrix.hpp

Go to the documentation of this file.

/*
Copyright (c) 2014-2016, Scott Zuyderduyn
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

The views and conclusions contained in the software and documentation are those
of the authors and should not be interpreted as representing official policies,
either expressed or implied, of the FreeBSD Project.
*/

//----------------------------------------------------------------------------
// matrix.hpp
// An STL-like structure that resides in video memory.
//
// Author: Scott D. Zuyderduyn, Ph.D. (scott.zuyderduyn@utoronto.ca)
//----------------------------------------------------------------------------

#pragma once
#ifndef ECUDA_MATRIX_HPP
#define ECUDA_MATRIX_HPP

#include <vector>

#include "global.hpp"
#include "algorithm.hpp"
#include "allocators.hpp"
#include "apiwrappers.hpp"
#include "memory.hpp"
#include "type_traits.hpp"

#include "model/device_contiguous_row_matrix.hpp"

namespace ecuda {

namespace impl {

template<typename T,class Alloc> class matrix_kernel_argument; // forward declaration

} // namespace impl

template< typename T, class Alloc=device_pitch_allocator<T>, class P=shared_ptr<T> >
class matrix : private model::device_contiguous_row_matrix< T, /*padded_ptr< T,*/P/* >*/ >
{

private:
    typedef model::device_contiguous_row_matrix< T, /*padded_ptr< T,*/P/* >*/ > base_type;

public:
    typedef typename base_type::value_type      value_type;      
    typedef Alloc                               allocator_type;  
    typedef typename base_type::size_type       size_type;       
    typedef typename base_type::difference_type difference_type; 
    typedef typename base_type::reference       reference;       
    typedef typename base_type::const_reference const_reference; 
    typedef typename base_type::pointer         pointer;         
    typedef typename make_const<pointer>::type  const_pointer;   

    typedef typename base_type::row_type          row_type;          
    typedef typename base_type::column_type       column_type;       
    typedef typename base_type::const_row_type    const_row_type;    
    typedef typename base_type::const_column_type const_column_type; 

    typedef typename base_type::iterator               iterator;               
    typedef typename base_type::const_iterator         const_iterator;         
    typedef typename base_type::reverse_iterator       reverse_iterator;       
    typedef typename base_type::const_reverse_iterator const_reverse_iterator; 

    typedef       impl::matrix_kernel_argument<T,Alloc> kernel_argument;       
    typedef const impl::matrix_kernel_argument<T,Alloc> const_kernel_argument; 

private:
    allocator_type allocator;
    //template<typename U,class Alloc2> class device_matrix;
    template<typename U,class Alloc2,class Q> friend class matrix;

protected:

    template<typename U>
    __HOST__ __DEVICE__ matrix( const matrix<T,Alloc,U>& src, ecuda::true_type ) : base_type( unmanaged_cast(src.get_pointer()), src.number_rows(), src.number_columns() ), allocator(src.allocator) {}

    template<typename U>
    __HOST__ __DEVICE__ matrix& shallow_assign( const matrix<T,Alloc,U>& other )
    {
        base_type::get_pointer() = unmanaged_cast(other.get_pointer());
        allocator = other.allocator;
        return *this;
    }

private:
    __HOST__ void init()
    {
        if( number_rows() && number_columns() ) {
            // TODO: this is unfortunate - have to get a padded_ptr from the allocator, unwrap it and
            //       give it to shared_ptr, and then rewrap it in a padded_ptr with the same attributes
            //       as the original - the device_contiguous_row_matrix second template parameter which
            //       enforces a padded_ptr of some type is the reason
            typename Alloc::pointer p = get_allocator().allocate( number_columns(), number_rows() );
            shared_ptr<value_type> sp( naked_cast<typename ecuda::add_pointer<value_type>::type>(p) );
            padded_ptr< value_type, shared_ptr<value_type> > pp( sp, p.get_pitch() ); //, p.get_width(), sp );
            base_type base( pp, number_rows(), number_columns() );
            base_type::swap( base );
        }
    }

public:
    __HOST__ matrix( const size_type numberRows=0, const size_type numberColumns=0, const value_type& value = value_type(), const allocator_type& allocator = allocator_type() ) :
        base_type( pointer(), numberRows, numberColumns ),
        allocator(allocator)
    {
        init();
        if( size() ) ecuda::fill( begin(), end(), value );
    }

    __HOST__ matrix( const matrix& src ) :
        base_type( pointer(), src.number_rows(), src.number_columns() ),
        allocator(src.get_allocator())
// TODO: this is broken due to some complaints from stdlib about making an iterator to a void pointer
//      allocator(std::allocator_traits<allocator_type>::select_on_container_copy_construction(src.get_allocator()))
    {
        init();
        if( size() ) ecuda::copy( src.begin(), src.end(), begin() );
    }

    __HOST__ matrix( const matrix& src, const allocator_type& alloc ) :
        base_type( pointer(), src.number_rows(), src.number_columns() ),
        allocator(alloc)
    {
        init();
        if( size() ) ecuda::copy( src.begin(), src.end(), begin() );
    }

    __HOST__ matrix& operator=( const matrix& other )
    {
        allocator = other.allocator;
        if( number_rows() != other.number_rows() || number_columns() != other.number_columns() )
            resize( other.number_rows(), other.number_columns() );
        if( size() ) ecuda::copy( other.begin(), other.end(), begin() );
        return *this;
    }

    #ifdef ECUDA_CPP11_AVAILABLE
    __HOST__ matrix( matrix&& src ) : base_type() { swap(src); }

    __HOST__ matrix& operator=( matrix&& src )
    {
        swap(src);
        return *this;
    }
    #endif

    __HOST__ __DEVICE__ inline iterator begin() __NOEXCEPT__ { return base_type::begin(); }

    __HOST__ __DEVICE__ inline iterator end() __NOEXCEPT__ { return base_type::end(); }

    __HOST__ __DEVICE__ inline const_iterator begin() const __NOEXCEPT__ { return base_type::begin(); }

    __HOST__ __DEVICE__ inline const_iterator end() const __NOEXCEPT__ { return base_type::end(); }

    __HOST__ __DEVICE__ inline reverse_iterator rbegin() __NOEXCEPT__ { return base_type::rbegin(); }

    __HOST__ __DEVICE__ inline reverse_iterator rend() __NOEXCEPT__ { return base_type::rend(); }

    __HOST__ __DEVICE__ inline const_reverse_iterator rbegin() const __NOEXCEPT__ { return base_type::rbegin(); }

    __HOST__ __DEVICE__ inline const_reverse_iterator rend() const __NOEXCEPT__ { return base_type::rend(); }

    #ifdef ECUDA_CPP11_AVAILABLE
    __HOST__ __DEVICE__ inline const_iterator         cbegin()  const __NOEXCEPT__ { return base_type::cbegin();  }
    __HOST__ __DEVICE__ inline const_iterator         cend()    const __NOEXCEPT__ { return base_type::cend();    }
    __HOST__ __DEVICE__ inline const_reverse_iterator crbegin()       __NOEXCEPT__ { return base_type::crbegin(); }
    __HOST__ __DEVICE__ inline const_reverse_iterator crend()         __NOEXCEPT__ { return base_type::crend();   }
    #endif

    __HOST__ __DEVICE__ inline size_type size() const __NOEXCEPT__ { return base_type::size(); }

    __HOST__ __DEVICE__ __CONSTEXPR__ inline size_type max_size() const __NOEXCEPT__ { return std::numeric_limits<size_type>::max(); }

    __HOST__ __DEVICE__ inline size_type number_rows() const __NOEXCEPT__ { return base_type::number_rows(); }

    __HOST__ __DEVICE__ inline size_type number_columns() const __NOEXCEPT__ { return base_type::number_columns(); }

    __HOST__ void resize( const size_type newNumberRows, const size_type newNumberColumns, const value_type& value = value_type() )
    {
        if( number_rows() == newNumberRows && number_columns() == newNumberColumns ) return; // no resize needed
        // create new model
        matrix newMatrix( newNumberRows, newNumberColumns, value, get_allocator() );
        for( size_type i = 0; i < std::min(number_rows(),newNumberRows); ++i ) {
            const_row_type oldRow = get_row(i);
            row_type newRow = newMatrix.get_row(i);
            ecuda::copy( oldRow.begin(), oldRow.begin()+std::min(number_columns(),newNumberColumns), newRow.begin() );
        }
        swap( newMatrix );
    }

    __HOST__ __DEVICE__ inline bool empty() const __NOEXCEPT__ { return !number_rows() || !number_columns(); }

    __HOST__ __DEVICE__ inline row_type get_row( const size_type rowIndex ) { return base_type::get_row(rowIndex); }

    __HOST__ __DEVICE__ inline const_row_type get_row( const size_type rowIndex ) const { return base_type::get_row(rowIndex); }

    __HOST__ __DEVICE__ inline column_type get_column( const size_type columnIndex ) { return base_type::get_column(columnIndex); }

    __HOST__ __DEVICE__ inline const_column_type get_column( const size_type columnIndex ) const { return base_type::get_column(columnIndex); }

    __DEVICE__ inline reference at( size_type rowIndex, size_type columnIndex )
    {
        if( rowIndex >= number_rows() || columnIndex >= number_columns() ) {
            #ifndef __CUDACC__
            throw std::out_of_range( EXCEPTION_MSG("ecuda::matrix::at() row and/or column index parameter is out of range") );
            #else
            // this strategy is taken from:
            // http://stackoverflow.com/questions/12521721/crashing-a-kernel-gracefully
            ecuda::threadfence();
            //__threadfence();
            asm("trap;");
            #endif
        }
        return base_type::at( rowIndex, columnIndex );
    }

    __DEVICE__ inline const_reference at( size_type rowIndex, size_type columnIndex ) const
    {
        if( rowIndex >= number_rows() || columnIndex >= number_columns() ) {
            #ifndef __CUDACC__
            throw std::out_of_range( EXCEPTION_MSG("ecuda::matrix::at() row and/or column index parameter is out of range") );
            #else
            // this strategy is taken from:
            // http://stackoverflow.com/questions/12521721/crashing-a-kernel-gracefully
            ecuda::threadfence();
            //__threadfence();
            asm("trap;");
            #endif
        }
        return base_type::at( rowIndex, columnIndex );
    }

    __DEVICE__ inline reference operator()( const size_type rowIndex, const size_type columnIndex ) { return base_type::at(rowIndex,columnIndex); }

    __DEVICE__ inline const_reference operator()( const size_type rowIndex, const size_type columnIndex ) const { return base_type::at(rowIndex,columnIndex); }

    __HOST__ __DEVICE__ inline row_type operator[]( const size_type rowIndex ) { return get_row(rowIndex); }

    __HOST__ __DEVICE__ inline const_row_type operator[]( const size_type rowIndex ) const { return get_row(rowIndex); }

    __DEVICE__ inline reference front() { return base_type::at(0,0); }

    __DEVICE__ inline reference back() { return base_type::at(number_rows()-1,number_columns()-1); }

    __DEVICE__ inline const_reference front() const { return base_type::at(0,0); }

    __DEVICE__ inline const_reference back() const { return base_type::at(number_rows()-1,number_columns()-1); }

    __HOST__ __DEVICE__ inline pointer data() __NOEXCEPT__ { return base_type::get_pointer(); }

    __HOST__ __DEVICE__ inline const_pointer data() const __NOEXCEPT__ { return base_type::get_pointer(); }

    __HOST__ __DEVICE__ inline void fill( const value_type& value ) { if( !empty() ) ecuda::fill( begin(), end(), value ); }

    __HOST__ __DEVICE__ inline void swap( matrix& other ) { base_type::swap( other ); }

    __HOST__ inline allocator_type get_allocator() const { return allocator; }

    template<class Alloc2>
    __HOST__ __DEVICE__ inline bool operator==( const matrix<value_type,Alloc2>& other ) const
    {
        return ecuda::equal( begin(), end(), other.begin() );
    }

    template<class Alloc2>
    __HOST__ __DEVICE__ inline bool operator!=( const matrix<value_type,Alloc2>& other ) const { return !operator==(other); }

    template<class Alloc2>
    __HOST__ __DEVICE__ inline bool operator<( const matrix<value_type,Alloc2>& other ) const
    {
        return ecuda::lexicographical_compare( begin(), end(), other.begin(), other.end() );
    }

    template<class Alloc2>
    __HOST__ __DEVICE__ inline bool operator>( const matrix<value_type,Alloc2>& other ) const
    {
        return ecuda::lexicographical_compare( other.begin(), other.end(), begin(), end() );
    }

    template<class Alloc2>
    __HOST__ __DEVICE__ inline bool operator<=( const matrix<value_type,Alloc2>& other ) const { return !operator>(other); }

    template<class Alloc2>
    __HOST__ __DEVICE__ inline bool operator>=( const matrix<value_type,Alloc2>& other ) const { return !operator<(other); }

};

namespace impl {

template< typename T, class Alloc> //=device_pitch_allocator<T> >
class matrix_kernel_argument : public matrix<T,Alloc,typename ecuda::add_pointer<T>::type>
{

private:
    typedef matrix<T,Alloc,typename ecuda::add_pointer<T>::type> base_type;

public:
    template<class P>
    __HOST__ matrix_kernel_argument( const matrix<T,Alloc,P>& src ) : base_type( src, ecuda::true_type() ) {}

    __HOST__ __DEVICE__ matrix_kernel_argument( const matrix_kernel_argument& src ) : base_type( src, ecuda::true_type() ) {}

    template<class P>
    __HOST__ matrix_kernel_argument& operator=( const matrix<T,Alloc,P>& src )
    {
        base_type::shallow_assign( src );
        return *this;
    }

    #ifdef ECUDA_CPP11_AVAILABLE
    matrix_kernel_argument( matrix_kernel_argument&& src ) : base_type(std::move(src)) {}

    matrix_kernel_argument& operator=( matrix_kernel_argument&& src )
    {
        base_type::operator=(std::move(src));
        return *this;
    }
    #endif

};

} // namespace impl

template<typename T,class Alloc1,class Alloc2>
__HOST__ void matrix_copy( matrix<T,Alloc1>& dest, const matrix<T,Alloc2>& src, typename matrix<T,Alloc2>::size_type offsetRow=0, typename matrix<T,Alloc2>::size_type offsetColumn=0 ) {
    typedef typename matrix<T,Alloc2>::size_type size_type;
    const size_type nr = std::min( dest.number_rows()   , src.number_rows()-offsetRow       );
    const size_type nc = std::min( dest.number_columns(), src.number_columns()-offsetColumn );
    for( size_type i = 0; i < nr; ++i ) {
        typename matrix<T,Alloc1>::row_type destRow = dest[i];
        typename matrix<T,Alloc2>::const_row_type srcRow = src[i+offsetRow];
        ::ecuda::copy( src.begin(), src.end(), dest.begin() );
    }
}

template<typename T,class Alloc1,class Alloc2>
__HOST__ void matrix_swap(
    matrix<T,Alloc1>& mat1,
    matrix<T,Alloc2>& mat2,
    typename matrix<T,Alloc1>::size_type numberRows=0, typename matrix<T,Alloc1>::size_type numberColumns=0,
    typename matrix<T,Alloc1>::size_type offsetRow1=0, typename matrix<T,Alloc1>::size_type offsetColumn1=0,
    typename matrix<T,Alloc2>::size_type offsetRow2=0, typename matrix<T,Alloc2>::size_type offsetColumn2=0
)
{
    if( (offsetRow1+numberRows) > mat1.number_rows() ) throw std::out_of_range( EXCEPTION_MSG("ecuda::matrix_swap() specified row subset of mat1 is out of bounds") );
    if( (offsetRow2+numberRows) > mat2.number_rows() ) throw std::out_of_range( EXCEPTION_MSG("ecuda::matrix_swap() specified row subset of mat2 is out of bounds" ) );
    if( (offsetColumn1+numberColumns) > mat1.number_columns() ) throw std::out_of_range( EXCEPTION_MSG("ecuda::matrix_swap() specified column subset of mat1 is out of bounds") );
    if( (offsetColumn2+numberColumns) > mat2.number_columns() ) throw std::out_of_range( EXCEPTION_MSG("ecuda::matrix_swap() specified column subset of mat2 is out of bounds") );
    std::vector< T, host_allocator<T> > stagingMemory( numberColumns );
    typedef typename matrix<T,Alloc1>::size_type size_type;
    for( size_type i = 0; i < numberRows; ++i ) {
        typename matrix<T,Alloc1>::row_type row1 = mat1[offsetRow1+i];
        typename matrix<T,Alloc2>::row_type row2 = mat1[offsetRow2+i];
        stagingMemory.assign( row1.begin()+offsetColumn1, row1.begin()+(offsetColumn1+numberColumns) );
        ecuda::copy( row2.begin()+offsetColumn2, row2.begin()+(offsetColumn2+numberColumns), row1.begin()+offsetColumn1 );
        ecuda::copy( stagingMemory.begin(), stagingMemory.end(), row2.begin()+offsetColumn2 );
    }
}

template<typename T,class Alloc>
__HOST__ void matrix_transpose(
    matrix<T,Alloc>& src
)
{
    if( src.empty() ) return;
    std::vector< T, host_allocator<T> > stagingMemory( src.number_columns() ); // stage a single row
    std::vector<T> hostMatrix( src.size() );
    for( typename matrix<T,Alloc>::size_type i = 0; i < src.number_rows(); ++i ) {
        ecuda::copy( src[i].begin(), src[i].end(), stagingMemory.begin() ); // copy row
        typename std::vector< T, host_allocator<T> >::const_iterator srcElement = stagingMemory.begin();
        for( std::size_t j = 0; j < src.number_columns(); ++j, ++srcElement ) hostMatrix[j*src.number_rows()+i] = *srcElement; // transpose
    }
    src.resize( src.number_columns(), src.number_rows() ); // resize destination matrix
    ecuda::copy( hostMatrix.begin(), hostMatrix.end(), src.begin() );
    //typename std::vector<T>::const_iterator srcRow = hostMatrix.begin();
    //for( typename matrix<T,Alloc>::size_type i = 0; i < src.number_rows(); ++i, srcRow += src.number_columns() ) {
    //  ecuda::copy( srcRow, srcRow+src.number_columns(), src[i].begin() );
    //}
}

} // namespace ecuda

#endif