Re: [ublas] blocked products and relation UBLAS/GLAS

From: Gunter Winkler (guwi17_at_[hidden])
Date: 2006-02-16 04:53:46

Riccardo Rossi wrote:

>Dear Gunther,
>     first of all thanks for your quick answer, i have however a doubt 
>concerning your proposal:
>  
>
I will comment on your mail later. But maybe you can already have a look 
at some old code of mine. I don't think it compiles any more - but I 
think about it later today.

mfg
Gunter


#include <iostream>

#define BOOST_UBLAS_USE_ET

#include <boost/timer.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>

using std::size_t;
using std::cout;
using std::endl;

// define matrix/vector types

typedef boost::numeric::ublas::bounded_matrix<double, 3, 3>  TINY;
typedef boost::numeric::ublas::matrix<TINY>                  BIG;

typedef boost::numeric::ublas::bounded_vector<double, 3>  TINY_VEC;
typedef boost::numeric::ublas::vector<TINY_VEC>           BIG_VEC;

// define new promote traits

namespace boost { namespace numeric { namespace ublas {

    template<>
    struct promote_traits<TINY, TINY_VEC> {
        typedef TINY_VEC promote_type;
    };

    template<>
    struct promote_traits<TINY, double> {
        typedef TINY promote_type;
    };

    template<>
    struct promote_traits<double, TINY> {
        typedef TINY promote_type;
    };

    BOOST_UBLAS_INLINE
        TINY_VEC
        operator * (const TINY& A, const TINY_VEC& x)
        {
          TINY_VEC tmp;
          tmp = prod(A,x);
          return tmp;
        };

  /* helper for matrix_unary2::operator()(...)
        // Element access
        BOOST_UBLAS_INLINE
        const_reference operator () (size_type i, size_type j) const {
            return functor_type () (trans(e_ (j, i))); 
        }
   */  
  inline
  double trans(const double x) 
  { 
        return x; 
  }

}}}

using namespace boost::numeric::ublas;

int main(size_t argc, char *argv[])
{
  size_t size = 20;

  if (argc > 1)
        size = ::atoi (argv [1]);

  TINY    tiny;
  BIG     big(5,5);

  tiny(0,0) = 1;
  tiny(0,1) = 2;
  tiny(0,2) = 3;
  tiny(1,0) = 4;
  tiny(2,0) = 5;
  
  cout << "tiny: " << tiny << endl;

  big(0,0) = tiny;
  big(1,1) = tiny;
  big(2,2) = tiny;
  big(3,4) = tiny;
  big(4,3) = tiny;

  cout << "big: " << big << endl;
  
  TINY_VEC  tv;

  tv(0) = 1.0;
  tv(1) = 1.0;
  tv(2) = 1.0;

  BIG_VEC   bv(scalar_vector<TINY_VEC>(5,tv));

  cout << "tiny_vector: " << tv << endl;
  cout << "big_vector: " << bv << endl;

  cout << "tiny + tiny = " << (tiny + tiny) << endl;
  cout << "big + big   = " << (big + big) << endl;

  cout << "tiny - tiny = " << (tiny - tiny) << endl;
  cout << "big - big   = " << (big - big) << endl;

  cout << "2.0 * tiny = " << (2.0 * tiny) << endl;
  cout << "2.0 * big  = " << (2.0 * big) << endl;

  cout << "tiny * tiny_vector = " << prod(tiny, tv) << endl;
  cout << "big * big_vector   = " << prod(big, bv) << endl;

  cout << "trans(tiny) = " << trans(tiny) << endl;
  cout << "trans(big)  = " << trans(big) << endl;

  return 0;

}