|
|
Ublas : |
Subject: Re: [ublas] matrix_range for matrix
From: oswin krause (oswin.krause_at_[hidden])
Date: 2013-07-06 13:35:53
Hi,
of course the overloaded operator= and missing accessors of matrix_range
are a nuisance. But this is nothing one could not work-around with a few
lines of code. Be creative :). For example interface_ does not need to
be a member, it could also be a temporary object. Or we could just ask
for a new function so that we can do at least something like
interface_.copy_proxy(subrange(matrix,0,size1,0,size2));
Greetings,
Oswin
On 06.07.2013 16:54, Petros wrote:
> Hi Oswin,
> The problem with this construction, is that you cannot resize this
> matrix.
> Petros
>
> -----Original Message----- From: oswin krause Sent: Saturday, July 06,
> 2013 4:28 AM To: ublas_at_[hidden] Subject: Re: [ublas]
> matrix_range for matrix
> Hi,
>
> there is no need to dig into the internals of ublas. In your case you
> could implment your class in terms of ublas::matrix and
> ublas::matrix_range and than just forward the public interface of the
> range. This could look similar to this (just the idea, no compilable
> code...):
>
> class MatrixGeneralT{
> public:
> MatrixGeneralT(
> const size_t nbrRows,
> const size_t nbrCols = 0
> )
> : matrix_( determineLDA( nbrRows ), nbrCols )
> , interface_(subrange(matrix_,0,nbrRows,0,nbrCols))//connect
> the matrix with the subrange
> {}
> //...
> typedef typename matrix_range_type::iterator1 iterator1;
> typedef typename matrix_range_type::iterator2 iterator2;
> //....
> const_reference operator()(size_type i, size_type j)const{
> return interface_(i,j);
> }
> reference operator()(size_type i, size_type j){
> return interface_(i,j);
> }
> size_type size1()const{
> return interface_.size1();
> }
> size_type size2()const{
> return interface_.size1();
> }
>
> //rest in the same manner: begin1/begin2 end1/end2 find1/find2 and
> that should be everything needed i guess...
>
>
> private:
> internal_matrix_type matrix_;
> matrix_range_type interface_;
> };
>
> On 06.07.2013 00:34, Petros wrote:
>> Naso,
>>
>> Other than a)the extra constructor with the default nbrCols argument
>> and the b)namespace alias vs.
>> the keyword using ( actually the latter is much more loose and better
>> be avoided, since other things get
>> defined as well ) the two pieces of code are identical.
>> ( the c/tor uses size2_ before it is initialized - this is the
>> problem in icl as well, only icl does not
>> make a fuss about it if you don't use it ;-) ).
>>
>> I am using the intel 12.1 compiler with msvc10 on windows and c++0x.
>> I don't understand why gcc would have any issue (maybe forgot the
>> std0x flag?) w/it but I do
>> not build with gcc, anyway..
>>
>> As far as rethinking the problem, I have thought about it a lot and
>> about a year ago as well and
>> cannot see any other way, without getting into the ublas internals -
>> in the end of the day, I
>> am only a user of ublas.
>> The way I see it, I need my matrix for 2 things: its own public
>> interface - I can take care of it, and its
>> connectivity with the ET of ublas. The latter is in the form of
>> either assignments ( again I can take
>> care of it ) and how it behaves on the lhs of expressions. The prod
>> example, shows that the trick of the
>> auto cast probably does the job.
>>
>> I know it is not ideal, since it is uncontrollable and could
>> "collide" with something else. I am hopeful because, in reality,
>> it is no more than another matrix expression . I give it the look and
>> feel of a matrix_range but "package" together the
>> matrix it refers to.
>> Maybe somebody else, who will have a better idea, given the
>> constraints, will come up and present it
>> in the forum, or even better add it to the code. Until then ..
>>
>> Thank you for all your help,
>> Petros
>>
>>
>> -----Original Message----- From: Nasos Iliopoulos
>> Sent: Friday, July 05, 2013 4:44 PM
>> To: ublas_at_[hidden]
>> Subject: Re: [ublas] matrix_range for matrix
>>
>> Petro,
>> here is a compilable and working version. GCC 4.7. had many issues with
>> the previous code. What compiler are you using? You can probably carve
>> this to fit your needs but I would suggest to rethink your approach.
>>
>> Let me know if it works.
>>
>> Best,
>> -Nasos
>>
>> #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>
>>
>>
>> using namespace boost::numeric;
>>
>> /* 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 =
>> ldiv( 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
>> MatrixGeneralT self_type ;
>> typedef typename
>> ublas::matrix< _T, _L, _A > M ;
>> typedef
>> M internal_matrix_type ;
>> typedef
>> 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 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
>> const_matrix_range_type const_closure_type;
>> typedef
>> 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 :
>> size_type size1_ ;
>> size_type size2_ ;
>> internal_matrix_type matrix_ ;
>>
>> 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 )
>> , matrix_( determineLDA( nbrRows ), nbrCols )
>> {}
>> MatrixGeneralT(
>> const size_t nbrRows
>> )
>> :size1_( nbrRows )
>> ,size2_( nbrRows )
>> , matrix_( determineLDA( nbrRows ), nbrRows )
>> {}
>>
>> 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;
>> }
>>
>>
>> On 07/05/2013 03:40 PM, Petros wrote:
>>> 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]
>>>
>>> _______________________________________________
>>> 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: Oswin.Krause_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: Oswin.Krause_at_[hidden]