|
|
Ublas : |
Subject: Re: [ublas] matrix_range for matrix
From: Petros (pmamales_at_[hidden])
Date: 2013-07-05 14:36:41
Naso,
As you can see there is more than just an aligned allocator . Every matrix
column has to be aligned,
making the ordinary unbounded_array inadequate in its current form as
(general) matrix storage(unless the
distinction between leading dimension and size1 for fortran matrices is
made, in which case the matrix code will have
to get modified - in reality the indexing happens in the layout types - via
static functions, that are filled with the array sizes,
obviously to accommodate matrix_container's and matrix proxies at once ).
( Of course allocators are always what they are and unbounded array is a
"glorified" c-array buffer, but I think you
know what I mean..)
Per your second point, since my problem is compilation, your request for
compilable code comes a bit tricky ;-))
However,
the good news is that I managed to make the matrix product compile with some
modifications on the const_closure_type
(btw, is there anywhere some documentation for these typedef's ? I
understand the unlikeness, but since these correspond
to design concepts ..) and the subsequent definition of iterators.
I will sent to the mailing list the basic ingredients of this kind of matrix
implementation shortly, in case anyone would be
interested, as a thank you to the list.
The assignment operator, through the auto cast, will not work however! I
suspect that c++ forbids it! The reason is that,
this way you would be able to bypass the rule that the assignment operator
is not definable outside the class scope !
Now, I can be totally wrong and some "guru" from the list w/c/should correct
me ;-). If there is a definite answer to this issue,
I would be obliged to anyone who would point it out to me.
For now, however, since the functions that have to be defined are only a
few, I can live with the present scaffolding.
Thank you very much for your help,
Petros P. Mamales
-----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]