#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))); 
        }
     */  

  template<class T>
  struct scalar_transpose:
    public scalar_unary_functor<T> {
    typedef const T & argument_type;
    typedef typename matrix_unary2_traits<const typename T::matrix_type, scalar_identity<typename T::value_type> >::expression_type result_type;

    static
    result_type apply (argument_type t) {
      return trans(t);
    }
  };

  // I hab to use 'big_trans' instead of 'trans' because the template matching
  // algorithm chose the wrong specialization/overload
  BOOST_UBLAS_INLINE
  matrix_unary2_traits<const BIG, scalar_transpose<BIG::value_type> >::result_type
  big_trans (const matrix_expression<BIG> &e) {
    typedef matrix_unary2_traits<const BIG, scalar_transpose<BIG::value_type> >::expression_type expression_type;
    return expression_type (e ());
  }

}}}

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.clear();
  tiny(0,0) = 1;
  tiny(0,1) = 2;
  tiny(0,2) = 3;
  tiny(1,1) = 4;
  tiny(2,1) = 5;
  
  cout << "tiny: " << tiny << endl;

  // WARINING: clear only works if the _default_ contructor of the value_type 
  // yields the zero element ( tiny() only creates _undefined_ values )
  big.clear();

  TINY tiny_zero(zero_matrix<double>(3,3));
  big = scalar_matrix<TINY>(big.size1(), big.size2(), tiny_zero);

  big(0,0) = tiny;
  big(1,0) = tiny;
  big(2,0) = tiny;
  big(3,0) = tiny;
  big(4,0) = 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 << "trans(tiny)    = " << trans(tiny) << endl;

  // this fails because trans() only transposes the big matrix
  // but not the single elements
  // read: they are modified by the scalar_identity functor 
  // see functional.hpp
  cout << "trans(big)     = " << trans(big) << endl;

  cout << "big_trans(big) = " << big_trans(big) << endl;

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

  // this product fails because the matrix product constructs
  // a temporary TINY_VEC using the "zero constructor"
  // (functional.hpp:852) TINY_VEC t(0);
  // which creates a bounded vector of size zero instead of
  // the expected 3-vector with value zero.
  cout << "big * big_vector   = " << prod(big, bv) << endl;

  return 0;

}