|
|
Ublas : |
Subject: Re: [ublas] matrix_range for matrix
From: Nasos Iliopoulos (nasos_i_at_[hidden])
Date: 2013-07-05 13:15:25
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]