|
|
Ublas : |
Subject: Re: [ublas] matrix_range for matrix
From: pmamales_at_[hidden]
Date: 2013-07-05 13:04:22
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]