Re: matrix_range for matrix
Nasos Iliopoulos <[email protected]>
| Newsgroups | gmane.comp.lib.boost.ublas |
|---|---|
| Message-ID | <[email protected]> |
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, [email protected] 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 <[email protected]> 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, [email protected] 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 >>> [email protected] >>> http://lists.boost.org/mailman/listinfo.cgi/ublas >>> Sent to: [email protected] >> _______________________________________________ >> ublas mailing list >> [email protected] >> http://lists.boost.org/mailman/listinfo.cgi/ublas >> Sent to: [email protected] > _______________________________________________ > ublas mailing list > [email protected] > http://lists.boost.org/mailman/listinfo.cgi/ublas > Sent to: [email protected] _______________________________________________ ublas mailing list [email protected] http://lists.boost.org/mailman/listinfo.cgi/ublas Sent to: [email protected]