|
|
Ublas : |
Subject: [ublas] matrix_range for matrix
From: pmamales_at_[hidden]
Date: 2013-07-05 10:51:02
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;
}