|
|
Ublas : |
Subject: Re: [ublas] matrix_range for matrix
From: Petros (pmamales_at_[hidden])
Date: 2013-07-05 15:40:08
Naso,
This is the skeleton of the code, in case someone needs a starting point.
Please, rest assured that there are shortcomings and pitfalls but vary
likely
a better starting point than nothing:
//boost
#include <boost/numeric/ublas/storage.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
//#include <boost/numeric/ublas/triangular.hpp>
namespace ublas = boost::numeric::ublas;
/* Really a+ boost general matrix_range to
a properly-lda-sized boost matrix it
holds internally.
*/
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 = true
) {
static const size_t MKL_LD_DIVSOR_IN_BYTES = 32 ;
static const size_t MKL_LD_DIVISOR_DISALLOWED = 2028 ;
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( MKL_LD_DIVSOR_IN_BYTES ) ) ;
//if not, resize lda so that this is true.
if ( qr.rem != 0 )
bytesNbrRows =
( qr.quot + 1 ) * MKL_LD_DIVSOR_IN_BYTES ;
// but avoid the bad dimension ( e.g. 2048)
while ( bytesNbrRows % MKL_LD_DIVISOR_DISALLOWED == 0 )
bytesNbrRows += MKL_LD_DIVSOR_IN_BYTES ;
return
( ld_ = bytesNbrRows / sizeof( _T ) ) ;
}
public:
typedef typename
MatrixGeneralT< _T, _L, _A> self_type ;
typedef typename
ublas::matrix< _T, _L, _A > M ;
typedef typename
M internal_matrix_type ;
typedef typename
const M const_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 M matrix_type ;
typedef typename M::size_type size_type ;
typedef typename M::difference_type difference_type ;
typedef typename M::value_type value_type ;
typedef typename M::const_reference const_reference ;
typedef typename std::conditional<
std::is_const<M>::value,
typename M::const_reference,
typename M::reference
>::type reference
;
typedef typename std::conditional<
std::is_const<M>::value ,
typename M::const_closure_type,
typename M::closure_type
>::type
matrix_closure_type ;
typedef ublas::basic_range<
size_type
, difference_type
> range_type ;
typedef typename
const_matrix_range_type
const_closure_type;
typedef typename
matrix_range_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;
typedef typename
matrix_range_type::iterator1 iterator1 ;
typedef typename
matrix_range_type::iterator2 iterator2 ;
typedef typename
const_matrix_range_type::iterator1
const_iterator1 ;
typedef typename
const_matrix_range_type::iterator2
const_iterator2 ;
private :
typename internal_matrix_type matrix_ ;
typename size_type size1_ ;
typename size_type size2_ ;
public :
MatrixGeneralT( bool optimalLDA = true )
:size1_(0)
,size2_(0)
, matrix_(0,0)
{}
MatrixGeneralT(
const size_t nbrRows,
const size_t nbrCols = 0
)
:size1_( nbrRows )
,size2_( nbrCols == 0 ? size1_ : nbrCols )
, matrix_( determineLDA( nbrRows ), nbrCols )
{}
virtual ~MatrixGeneralT()
{}
// auto casting
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_ )
) ;
}
// matrix_expression assignment :
template <typename _E >
self_type &
operator = ( typename ublas::matrix_expression<_E> const & me ) {
_E const & e = me() ;
matrix_.resize( determineLDA( e.size1(), true ), e.size2() ) ;
size1_ = e.size1() ;
size2_ = e.size2() ;
matrix_range_type
(
matrix_
, ublas::range( 0, size1_ )
, ublas::range( 0, size2_ )
)
=
me() ;
return
*this ;
}
value_type & operator() ( const size_t i, const size_t j ) {
return
matrix_.data()[ i + ld_ * j ] ;
}
value_type const & operator() ( const size_t i, const size_t j ) const {
return
matrix_.data()[ i + ld_ * j ] ;
}
} ;
int main() {
typedef MatrixGeneralT< double > matrix ;
matrix m(3,3) ;
ublas::vector< double > x( 3 ) ;
x[0] = 1; x[1] = 2 ; x[2] = 3;
m = ublas::outer_prod( x,x ) ;
matrix m1( 3, 4 ) ;
matrix m2( 4, 5 ) ;
matrix m3( 3, 5 ) ;
for ( size_t i = 0; i != 3 ; ++i )
for ( size_t j = 0; j != 4 ; ++j )
m1(i,j) = (i+1) + (j+1)*10 ;
for ( size_t i = 0; i != 4 ; ++i )
for ( size_t j = 0; j != 5 ; ++j )
m2(i,j) = (i+1) + (j+1)*10 ;
m3 = ublas::prod( m1, m2 ) ;
return
1;
}
-----Original Message-----
From: Nasos Iliopoulos
Sent: Friday, July 05, 2013 1:15 PM
To: ublas_at_[hidden]
Subject: Re: [ublas] matrix_range for matrix
Two things:
As part of GSoC Sathyam is building vector/Matrix view classes, along
with aligned allocators, that might help you get a much more clean
solution when interfacing with external code. His vector code and
aligned allocators are already in good shape and I expect we will get
them into uBlas soon (as in 3 months).
Additionally if you could provide some compilable code it would help
debug the problem.
Best,
-Nasos
On 07/05/2013 01:04 PM, pmamales_at_[hidden] wrote:
> Naso,
> Thx for getting back to me.
>
> A clarification : the matrix_range'ness of my matrix is coming into play
> only through the auto cast operator!
> It is really smoke and mirrors. Tried to "pipe" my class into ublas by the
> typedefs and by deriving from
> matrix_expression.
>
> The auto cast operation was also failing me for the assignment of a
> ublas::matrix_expression.
> So, I said, fine, I will "view" the general matrix (matrix_) as a range
> and perform the assignment
> there. What is returned is the object itself, not the "view".
> Now why the assignment was failing me, I do not know either ;-)). But
> assignments are more tricky
> and I thought of the ublas::prod test.
>
> All this is a hackery, unfortunately. A clean solution, i.e. one that
> would provide with a
> leading dimension, decoupling the (fortran-layout matrices) leading
> dimension from size1() is not present.
> Very likely, ublas designers wanted this behavior to be taken care by
> matrix_range, which is the more
> appropriate c++ way.
>
> However, LAPACK has these "oddities" to account for the FORTRAN
> limitations ( at least back then ).
> You see, I use ublas structures as containers, for lapack operations -
> speed is essential in what I do.
> Occasionally, and for all kinds of test calculations ( or even in the
> private area of my clases,
> for at least temporary purposes) I would like my struct's to allow for the
> ublas operations.
>
> I am stuck with this for a couple of days now, and try to find a clean way
> out.
> The auto cast operation seemed like a good candidate. However, I get
> problems that I did not anticipate
> and do not understand..
> If you or anyone else could point to a better "take" I "am all ears" ;-)
>
> Thank you for your help
> Petros
>
> ps: Very likely, everybody understands that, I forgot to include the ublas
> headers and the namespace alias,
> in the code snippet I sent. Apologies..
>
> ---- Nasos Iliopoulos <nasos_i_at_[hidden]> wrote:
>> Petro, Excuse me, but I hit send before finishing the e-mail:
>>
>> What are you 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 ;
>> }
>>
>> it looks that you are creating a temporary that dies before ever being
>> assigned to the encapsulated matrix_renge
>>
>> If this is the case than the matrix_range will be null, and hence trying
>> to write to it will probably crash everything.
>>
>> -Nasos
>>
>>
>> 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]
>> _______________________________________________
>> ublas mailing list
>> ublas_at_[hidden]
>> http://lists.boost.org/mailman/listinfo.cgi/ublas
>> Sent to: pmamales_at_[hidden]
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: athanasios.iliopoulos.ctr.gr_at_[hidden]
_______________________________________________
ublas mailing list
ublas_at_[hidden]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: pmamales_at_[hidden]