|
|
Ublas : |
Subject: Re: [ublas] matrix_range for matrix
From: Nasos Iliopoulos (nasos_i_at_[hidden])
Date: 2013-07-05 11:51:44
Petro,
I am not sure what you are trying to achieve here:
template <typename _me>
matrix_ge_self_type &
operator = (_me const & me )
{
matrix_range_type
(
matrix_
, ublas::range( 0, size1_ )
, ublas::range( 0, size2_ )
) = me ;
return
*this ;
}
On 07/05/2013 10:51 AM, pmamales_at_[hidden] wrote:
> In order to accommodate for special memory layout for FORTRAN-type matrices I need to use a "matrix" which is really a ublas::matrix_range on a
> ublas::matrix ( memory has to be aligned on 32-byte memory and actually every column of the matrix has to use this kind of boundary ).
> In order to do this, I am defining my own ublas-conformant type (called MatrixGeneralT and declared as a ublas::matrix_expression) which has the same template arguments as the ublas::matrix and defines
> the type-"interface" of a matrix range. To make things "happen" I create an automatic conversion from my class to the ublas::matrix_range.
> However, when I try, as a first test, to calculate the product of 2 matrices, the matrix_matrix_binary constructor fails to convert my class to matrix_range !
> Unless, I have done something that c++ forbids ( which I don't think, since tried the "trick" with toy classes ) there must be something in ublas preventing this from
> working.
> Any ideas/input will be greatly appreciated.
> TIA for your help,
> Petros
> ps: I use intel compiler 12.1 on windows 7 64bit plugged into msvc2010 ( hence the stl "flavor" ) and C++0x.
>
> I supply here a pseudo code for anyone interested in it :
> #include <boost/numeric/ublas/matrix_expression.hpp>
> template < typename _T, typename _L = ublas::column_major, typename _A = ublas::unbounded_array<_T> >
> class MatrixGeneralT
> :public ublas::matrix_expression< MatrixGeneralT<_T, _L, _A> >
> {
> size_t ld_ ;
>
> size_t
> determineLDA (
> const size_t nRows
> , const bool mklLeadingDimension
> ) {
> if ( ! mklLeadingDimension )
> return ( ld_ = nRows ) ;
> ld_ = nRows ;
> // the nbr of bytes per lda should be divisible by MKL::LDA_DIVISOR_IN_BYTES.
> size_t bytesNbrRows =
> nRows * sizeof( _T ) ;
> const ldiv_t qr =
> div( long(bytesNbrRows), long( 32/*MKL::LDA_DIVISOR_IN_BYTES*/ ) ) ;
> //if not, resize lda so that this is true.
> if ( qr.rem != 0 )
> bytesNbrRows =
> ( qr.quot + 1 ) * 32 /*MKL::LDA_DIVISOR_IN_BYTES */;
> // but avoid the bad dimension ( e.g. 2048)
> while ( bytesNbrRows % MKL::LDA_DIVISOR_IN_BYTES_FORBITTEN == 0 )
> bytesNbrRows += MKL::LDA_DIVISOR_IN_BYTES ;
> return
> ( ld_ = bytesNbrRows / sizeof( _T ) ) ;
> }
>
> size_t size1_ ;
> size_t size2_ ;
>
> public:
> typedef typename
> _T value_type ;
> typedef typename
> ublas::matrix< _T, _L, _A > internal_matrix_type ;
> typedef typename
> ublas::matrix_range<
> internal_matrix_type
> > matrix_range_type ;
> typedef typename
> ublas::matrix_range<
> const internal_matrix_type
> > const_matrix_range_type ;
>
> typedef typename internal_matrix_type M ;
> typedef typename M matrix_type;
> typedef typename M::size_type size_type;
> typedef typename M::difference_type difference_type;
> typedef typename M::const_reference const_reference;
> typedef typename boost::mpl::if_<boost::is_const<M>,
> typename M::const_reference,
> typename M::reference>::type reference;
> typedef typename boost::mpl::if_<boost::is_const<M>,
> typename M::const_closure_type,
> typename M::closure_type>::type matrix_closure_type;
> typedef typename ublas::basic_range<size_type, difference_type> range_type;
>
> typedef typename matrix_range_type self_type ;
> typedef const self_type const_closure_type;
> typedef self_type closure_type;
> typedef typename ublas::storage_restrict_traits<typename M::storage_category,
> ublas::dense_proxy_tag>::storage_category storage_category;
> typedef typename M::orientation_category orientation_category;
>
> #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
> typedef indexed_iterator1<matrix_range<matrix_type>,
> typename subiterator1_type::iterator_category> iterator1;
> typedef indexed_iterator2<matrix_range<matrix_type>,
> typename subiterator2_type::iterator_category> iterator2;
> typedef indexed_const_iterator1<matrix_range<matrix_type>,
> typename const_subiterator1_type::iterator_category> const_iterator1;
> typedef indexed_const_iterator2<matrix_range<matrix_type>,
> typename const_subiterator2_type::iterator_category> const_iterator2;
> #else
> class const_iterator1;
> class iterator1;
> class const_iterator2;
> class iterator2;
> #endif
> typedef ublas::reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
> typedef ublas::reverse_iterator_base1<iterator1> reverse_iterator1;
> typedef ublas::reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
> typedef ublas::reverse_iterator_base2<iterator2> reverse_iterator2;
> private :
> typename internal_matrix_type matrix_ ;
> public :
> typedef MatrixGeneralT< _T, _L, _A> matrix_ge_self_type ;
> MatrixGeneralT( bool optimalLDA = true )
> :size1_(0), size2_(0)
> {}
> MatrixGeneralT( const size_t nbrRows,
> const size_t nbrCols = 0,
> bool optimalLDA = true
> )
> : size1_(nbrRows), size2_(nbrCols )
> , matrix_( determineLDA( nbrRows, optimalLDA), nbrCols )
> {}
> virtual ~MatrixGeneralT()
> {}
>
>
> template <typename _me>
> matrix_ge_self_type &
> operator = (_me const & me )
> {
> matrix_range_type
> (
> matrix_
> , ublas::range( 0, size1_ )
> , ublas::range( 0, size2_ )
> ) = me ;
> return
> *this ;
> }
> operator matrix_range_type()
> {
> return
> matrix_range_type
> (
> matrix_
> , ublas::range( 0, size1_ )
> , ublas::range( 0, size2_ )
> )
> }
> //operator const_matrix_range_type () const {
> // return
> // const_matrix_range_type
> // (
> // matrix_
> // , ublas::range( 0, size1_ )
> // , ublas::range( 0, size2_ )
> // )
> //}
>
> } ;
>
> int main() {
>
> typedef MatrixGeneralT< double > matrix ;
>
> matrix m(3,3) ;
> ublas::vector< double > x( 3 ) ;
> x[0] = 1; x[1] = 2 ; x[2] = 3;
> // matrix::matrix_range_type mr(m) ;
> m = ( ublas::outer_prod( x,x ) ) ;
>
>
> matrix m1( 3, 5 );
> matrix m2( 5, 4 ) ;
> matrix m3( 3, 4 ) ;
>
> m3 = ublas::prod( m1, m2 ) ;
> return
> 1;
> }
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: athanasios.iliopoulos.ctr.gr_at_[hidden]