[boost] Anyone interested in a generic loop (and loop unrolling) library ?

From: Hui Li (huili80_at_[hidden])
Date: 2008-06-28 12:20:50

Loops, in particular loop-unrolling, can be made generic, easy to write, and
accessible to everybody. With the help of Boost.Lambda (and a few other
boost libraries), we can easily construct arbitrary loops that are unrolled
at compile-time.

I haven't implemented the full library yet, in fact I've implemented only a
few very preliminary functionalities (hardly flawless), it serves only to
illustrate the ideas. From what they can do, a possible Boost.Loop library
looks very promising. However I feel that I need much help on the design of
such a library. The implementation could be straightforward (with the help
of so many good libraries in Boost), once a good design is found. Hopefully
there are interested people here so I won't have to work on my own.  (If
there already exists such a library, please let me know.)

In the example, you can see how easy it is to construct arbitrary (limited
only by Boost.Lamda) loops that are generically unrolled at compile time.
With an optimizing compiler, it can be as efficient as hand-written codes.
Not only that, but also can we make (almost) any linear-algebra libraries
work together without difficulty.  In the example, we don't need any extra
coding to make a c_array, std::vector, or even a column_view of a
ublas::matrix to work together.

In the following I give the full example code to illustrate what it's
capable of and how to use it. The header file, loop.hpp, is attached at the
end of this post.
The example has been tested with apple-gcc4.0.1 .

Thanks,

Hui

/**********************************************************************************************************/

#include <vector>
#include <iostream>
#include <iterator>
#include <functional>

#include <boost/array.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>

#include "loop.hpp"

using namespace std;
using namespace boost;
using namespace boost::lambda;
using namespace boost::numeric;
using namespace boost::numeric::ublas;

using namespace loop;

int idt(int i)
{
    return i;
}

int main(int argc, char* argv[])
{
    const size_t N = 8, M = 3;

    boost::array<double,N> v;
    std::vector<double> w(N);
    matrix<double> m(N,M);

    // All loops are unrolled at compile time

    // indtrans acts like an array whose individual element equals its index
    index_transform<double(int)> indtrans(&idt);

    // equivalent to for(i=0;i<N;i+=2){ v[i]=w[i]=1; }
    loop_over< linear_range<0,N,2> >( _1 = _2 = 1, v, w );

    // equivalent to for(i=1;i<N;i+=2){ v[i]=i; }
    loop_over< linear_range<1,N,2> >( _1 = _2, v, indtrans );

    // equivalent to for(i=1;i<N;i+=2){ w[i]=i; }
    loop_over< linear_range<1,N,2> >( _1 = _2, w, indtrans );
//    loop_over< linear_range<1,N,2> >( _1 = _2 = _3, v, w, indtrans ); // I
wonder why this line doesn't compile with gcc4??
    m.clear(); // set m to a zero matrix

    // equivalent to for(i=1;i<N;++i){ column(m,2)[i]=1+i*2; }
    loop_over< linear_range<0,N,1> >( _1 = 1+_2*2, loop::ref(column(m,2)),
indtrans );

    // equivalent to for(i=1;i<2;++i){ v[i]*=-2; } followed by
for(i=N-3;i<N;++i){ v[i]*=-2; }
    loop_over< link_range<linear_range<0,2,1>,linear_range<N-3,N,1> > >( _1
*= -2, v);

    // vertically print out v and w, side by side
    loop_over< linear_range<0,N,1> >( cout << _1 << constant('\t') << _2 <<
constant('\n'), v, w ); cout<<endl;
    // print out m
    cout << m << endl << endl;

    double s = 0;
    // calculate the inner-product of w and column(m,2) (when the column is
viewed as a vector/array)
    loop_over< linear_range<0,N,1> >( _1 += _2*_3,
loop::ref(loop::scalar(s)), w, loop::ref(column(m,2)) );
    cout << s << endl << endl;

    return 0;
}

/**********************************************************************************************************/

Here is the loop.hpp

/**********************************************************************************************************
 *  loop.hpp
 *  Created by Hui Li on 6/26/08.
 ********************************************************************************************************
*/

#ifndef LOOP_H
#define LOOP_H

#include <boost/static_assert.hpp>
#include <boost/ref.hpp>
#include <boost/lambda/lambda.hpp>
#include <boost/function.hpp>
#include <boost/function_types/function_arity.hpp>
#include <boost/mpl/if.hpp>
#include <boost/mpl/bool.hpp>
#include <boost/mpl/assert.hpp>

namespace loop {

// the empty range
struct empty_range {};

// a range starting from Start, ending at End, with stride Stride
template < int Start, int End, int Stride = 1 >
struct linear_range
{
    BOOST_MPL_ASSERT(( boost::mpl::bool_<(Start*Stride<End*Stride)> ));
    static const int first  = Start;

    typedef typename boost::mpl::if_ <    boost::mpl::bool_<
((Start+Stride)*Stride < End*Stride) >,
                                        linear_range < Start+Stride, End,
Stride >,
                                        empty_range    >::type
next_range;
};

// link two arbitrary ranges
template < typename LeftRange, typename RightRange >
struct link_range
{
    static const int first = LeftRange::first;

    typedef typename boost::mpl::if_ <    typename boost::is_same < typename
LeftRange::next_range, empty_range >::type,
                                        RightRange,
                                        link_range < typename
LeftRange::next_range, RightRange >
                                     >::type    next_range;
};

// implementation of the unrolled loop over a range
template < typename LambdaExpr, typename Range >
struct loop_Impl
{
    template < typename V1 >
    void operator()(LambdaExpr& expr, V1& v1)
    {
        expr(v1[Range::first]);
        loop_Impl < LambdaExpr, typename Range::next_range >()(expr,v1);
    }

    template < typename V1, typename V2 >
    void operator()(LambdaExpr& expr, V1& v1, V2& v2)
    {
        expr(v1[Range::first],v2[Range::first]);
        loop_Impl < LambdaExpr, typename Range::next_range >()(expr,v1,v2);
    }

    template < typename V1, typename V2, typename V3 >
    void operator()(LambdaExpr& expr, V1& v1, V2& v2, V3& v3)
    {
        expr(v1[Range::first],v2[Range::first],v3[Range::first]);
        loop_Impl < LambdaExpr, typename Range::next_range
>()(expr,v1,v2,v3);
    }
};

//    looping over an empty range does nothing
template < typename LambdaExpr >
struct loop_Impl < LambdaExpr, empty_range >
{
    template < typename V1 >
    void operator()(LambdaExpr&, V1){}

    template < typename V1, typename V2 >
    void operator()(LambdaExpr&, V1&, V2&){}

    template < typename V1, typename V2, typename V3 >
    void operator()(LambdaExpr&, V1&, V2&, V3&){}
};

// helper functions
template < typename Range, typename LambdaExpr, typename V1 >
void loop_over(LambdaExpr& expr, V1& v1)
{
    loop_Impl<LambdaExpr,Range>()(expr,v1);
}

template < typename Range, typename LambdaExpr, typename V1, typename V2 >
void loop_over(LambdaExpr& expr, V1& v1, V2& v2)
{
    loop_Impl<LambdaExpr,Range>()(expr,v1,v2);
}

template < typename Range, typename LambdaExpr, typename V1, typename V2,
typename V3 >
void loop_over(LambdaExpr& expr, V1& v1, V2& v2, V3& v3)
{
    loop_Impl<LambdaExpr,Range>()(expr,v1,v2,v3);
}

// map a call to operator[] to a unary function all
template < typename Function >
struct index_transform
{
    BOOST_STATIC_ASSERT((
boost::function_types::function_arity<Function>::value == 1 ));

    typedef boost::function<Function>        func_type;

    typename func_type::result_type    operator[](typename
func_type::argument_type x)
    {
        return func(x);
    }

    template < typename FuncTp >
    index_transform(FuncTp f):func(f){}

private:
    func_type func;
};

// force removal of const-ness of a reference
// not so elegant but should be okay
template < typename T >
T& ref(const T& t)
{
    return const_cast<T&>(t);
}

// independent of whatever the index might be, it always returns the
reference to the same variable
template < typename T >
struct scalar_Impl
{
    T& operator[](int) { return ref; }
    scalar_Impl(T& t):ref(t){}
private:
    T& ref;
};

// scalar helper
template < typename T >
scalar_Impl<T> scalar(T& t)
{
    return scalar_Impl<T>(t);
}

} // namespace loop

#endif

/**********************************************************************************************************/