Re: matrix_range for matrix

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: [email protected]
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, [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] 


_______________________________________________
ublas mailing list
[email protected]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [email protected]