matrix.h

Go to the documentation of this file.

//matrix.h
//describes 'matrix' and 'vector' template classes and related stuff
//file version 0.0 20120622
//100% hand crafted
#ifndef TOOLS_MATRIX_INCLUDED
#define TOOLS_MATRIX_INCLUDED
#include <algorithm>
#include <cmath>
#include <stdexcept>
#pragma warning(push, 4)

namespace nabla {
#ifdef NDEBUG
#define ASSERT(cond) (void)(sizeof(cond?0:0))
#else
#include <cassert>
#define ASSERT(cond) assert(cond)
#endif


template<typename type>
class traits
{
public:
    static const type unit; 
    static const type zero; 
};

template<typename type> const type traits<type>::unit = static_cast<type>(1);
template<typename type> const type traits<type>::zero = static_cast<type>(0);

namespace tag {
struct storage {};

struct reference {};


template<size_t complexity>
struct expression {};
}//namespace tag

struct rectangular {};

struct symmetric {};

struct lower_triangular {};
typedef lower_triangular lower_tr;  

struct upper_triangular {};
typedef upper_triangular upper_tr;  


template<typename value_t, typename expr_t, typename tag>
struct matrix_expression
{};


struct constrained_index
{
    const size_t value;
    constrained_index(size_t i) : value(i) {}
    operator size_t() { return value; }

private:
    constrained_index &operator=(const constrained_index&);
};

namespace util {
template<typename tag>
struct tag_complexity
{ enum { value = 0 }; };
template<size_t c>
struct tag_complexity<tag::expression<c> >
{ enum { value = c }; };
struct unconstrained_first_non_zero
{
    typedef size_t first_type;
    template<typename type>
    static first_type first(const type&, size_t) { return 0; }
};
struct unconstrained_last_non_zero_in_row
{
    typedef size_t last_type;
    template<typename type>
    static last_type last(const type &m, size_t) { return m.cols() - 1; }
};
struct unconstrained_last_non_zero_in_column
{
    typedef size_t last_type;
    template<typename type>
    static last_type last(const type &m, size_t) { return m.rows() - 1; }
};
struct constrained_first_non_zero
{
    typedef constrained_index first_type;
    template<typename type>
    static first_type first(const type&, size_t i) { return i; }
};
struct constrained_last_non_zero
{
    typedef constrained_index last_type;
    template<typename type>
    static last_type last(const type &, size_t i) { return i; }
};
template<typename shape>
struct non_zero_in_row : unconstrained_first_non_zero, unconstrained_last_non_zero_in_row
{
    //typedef size_t first_type;
    //typedef size_t last_type;
    //template<typename type>
    //static first_type first(const type&, size_t) { return 0; }
    //template<typename type>
    //static last_type last(const type &m, size_t) { return m.cols() - 1; }
};
template<typename shape>
struct non_zero_in_column : unconstrained_first_non_zero, unconstrained_last_non_zero_in_column
{
    //typedef size_t first_type;
    //typedef size_t last_type;
    //template<typename type>
    //static first_type first(const type&, size_t) { return 0; }
    //template<typename type>
    //static last_type last(const type &m, size_t) { return m.rows() - 1; }
};
template<>
struct non_zero_in_row<lower_triangular> : unconstrained_first_non_zero, constrained_last_non_zero
{
    //typedef size_t first_type;
    //typedef constrained_index last_type;
    //template<typename type>
    //static first_type first(const type&, size_t) { return 0; }
    //template<typename type>
    //static last_type last(const type&, size_t row) { return row; }
};
template<>
struct non_zero_in_column<lower_triangular> : constrained_first_non_zero, unconstrained_last_non_zero_in_column
{
    //typedef constrained_index first_type;
    //typedef size_t last_type;
    //template<typename type>
    //static first_type first(const type&, size_t column) { return column; }
    //template<typename type>
    //static last_type last(const type &m, size_t) { return m.rows() - 1; }
};
template<>
struct non_zero_in_row<upper_triangular> : constrained_first_non_zero, unconstrained_last_non_zero_in_row
{
    //typedef constrained_index first_type;
    //typedef size_t last_type;
    //template<typename type>
    //static first_type first(const type&, size_t row) { return row; }
    //template<typename type>
    //static last_type last(const type &m, size_t) { return m.cols() - 1; }
};
template<>
struct non_zero_in_column<upper_triangular> : unconstrained_first_non_zero, constrained_last_non_zero
{
    //typedef size_t first_type;
    //typedef constrained_index last_type;
    //template<typename type>
    //static first_type first(const type&, size_t) { return 0; }
    //template<typename type>
    //static last_type last(const type &, size_t column) { return column; }
};
}//namespace util


template<typename value_t, typename expr_t, typename tag_t, typename shape>
struct shaped_expression :
    public matrix_expression<value_t,
                             expr_t,
                             tag::expression<util::tag_complexity<tag_t>::value>
                             >
{
    friend
        typename util::non_zero_in_row<shape>::first_type
        first_non_zero_in_row(const expr_t &m, size_t row)
    { return util::non_zero_in_row<shape>::first(m, row); }
    friend
        typename util::non_zero_in_row<shape>::last_type
        last_non_zero_in_row(const expr_t &m, size_t row)
    { return util::non_zero_in_row<shape>::last(m, row); }
    friend
        typename util::non_zero_in_column<shape>::first_type
        first_non_zero_in_column(const expr_t &m, size_t column)
    { return util::non_zero_in_column<shape>::first(m, column); }
    friend
        typename util::non_zero_in_column<shape>::last_type
        last_non_zero_in_column(const expr_t &m, size_t column)
    { return util::non_zero_in_column<shape>::last(m, column); }
};



template<typename expr_t> inline
util::non_zero_in_row<rectangular>::first_type
first_non_zero_in_row(const expr_t &m, size_t row)
{ return util::non_zero_in_row<rectangular>::first(m, row); }
template<typename expr_t> inline
util::non_zero_in_row<rectangular>::last_type
last_non_zero_in_row(const expr_t &m, size_t row)
{ return util::non_zero_in_row<rectangular>::last(m, row); }
template<typename expr_t> inline
util::non_zero_in_column<rectangular>::first_type
first_non_zero_in_column(const expr_t &m, size_t column)
{ return util::non_zero_in_column<rectangular>::first(m, column); }
template<typename expr_t> inline
util::non_zero_in_column<rectangular>::last_type
last_non_zero_in_column(const expr_t &m, size_t column)
{ return util::non_zero_in_column<rectangular>::last(m, column); }

template<typename expression_type> struct expr_traits;

template<typename value_t, typename expr_t, typename tag_t>
struct expr_traits<matrix_expression<value_t, expr_t, tag_t> >
{
    typedef value_t value_type; 
    typedef expr_t type;    
    typedef tag_t tag;  
    typedef rectangular shape;  
    enum {
        complexity = util::tag_complexity<tag>::value   
    };
};

template<typename value_t, typename expr_t, typename tag_t, typename shape_t>
struct expr_traits<shaped_expression<value_t, expr_t, tag_t, shape_t> >
{
    typedef value_t value_type; 
    typedef expr_t type;    
    typedef tag_t tag;  
    typedef shape_t shape;  
    enum {
        complexity = util::tag_complexity<tag>::value   
    };
};




template<typename value_t, typename expr_t, typename tag> inline
expr_t &get_matrix(matrix_expression<value_t, expr_t, tag> &a)
{ return *static_cast<expr_t*>(&a); }


template<typename value_t, typename expr_t, typename tag> inline
const expr_t &get_matrix(const matrix_expression<value_t, expr_t, tag> &a)
{ return *static_cast<const expr_t*>(&a); }

namespace orientation {
struct undefined {};

struct column {};

struct row {};
}//namespace orientation




template<typename expr_t>
orientation::undefined orientation_type(expr_t*);   //do not define

namespace util {    //all this allows to determine the orientation of a given vector expression
char id(orientation::undefined);
char (&id(orientation::column))[2];
char (&id(orientation::row))[3];
template<size_t> struct interpret_orientation;
template<> struct interpret_orientation<1> { typedef orientation::undefined type; };
template<> struct interpret_orientation<2> { typedef orientation::column type; };
template<> struct interpret_orientation<3> { typedef orientation::row type; };
template<typename expr_t>
struct vector_orientation
{
    typedef typename interpret_orientation<sizeof(id(orientation_type((expr_t*)42)))>::type type;
};
}//namespace util

//forward declaration
template<typename value_t, typename expr_t, typename tag,
         typename orientation = typename util::vector_orientation<expr_t>::type
         >
struct vector_expression;


template<typename value_t, typename expr_t, typename tag, typename orientation>
struct vector_expression {};


template<typename iterator_type>
struct base_row_iterator
{
    const size_t row;
    size_t col;
    base_row_iterator(size_t r, size_t c) : row(r), col(c) {}
    iterator_type &operator++() { ++col; return *static_cast<iterator_type*>(this); }
private:
    base_row_iterator &operator=(const base_row_iterator&);
};

template<typename t>
struct row_iterator
{
    struct type : public base_row_iterator<type>
    {
        t &ref;
        type(t &m, size_t r, size_t c) : base_row_iterator(r, c), ref(m) { ASSERT(row<ref.rows()); }
        typename t::reference operator*() { ASSERT(col<ref.cols()); return ref(row, col); }
    };
    static type begin(t &v, size_t r, size_t c) { return type(v, r, c); }
};

template<typename t>
struct row_const_iterator
{
    struct type : public base_row_iterator<type>
    {
        const t &ref;
        type(const t &m, size_t r, size_t c) : base_row_iterator(r, c), ref(m) { ASSERT(row<ref.rows()); }
        typename t::const_reference operator*() const { ASSERT(col<ref.cols()); return ref(row, col); }
    };
    static type begin(const t &v, size_t r, size_t c) { return type(v, r, c); }
};

template<typename iterator_type>
struct base_column_iterator
{
    const size_t col;
    size_t row;
    base_column_iterator(size_t r, size_t c) : row(r), col(c) {}
    iterator_type &operator++() { ++row; return *static_cast<iterator_type*>(this); }
private:
    base_column_iterator &operator=(const base_column_iterator&);
};

template<typename t>
struct column_iterator
{
    struct type : public base_column_iterator<type>
    {
        t &ref;
        type(type &m, size_t c, size_t r) : base_column_iterator(r, c), ref(m) { ASSERT(col<ref.cols()); }
        typename t::reference operator*() { ASSERT(row<ref.rows()); return ref(row, col); }
    };
    static type begin(const t &v, size_t c, size_t r) { return type(v, c, r); }
};

template<typename t>
struct column_const_iterator
{
    struct type : public base_column_iterator<type>
    {
        const t &ref;
        type(const t &m, size_t c, size_t r) : base_column_iterator(r, c), ref(m) { ASSERT(col<ref.cols()); }
        typename t::const_reference operator*() const { ASSERT(row<ref.rows()); return ref(row, col); }
    };
    static type begin(const t &v, size_t c, size_t r) { return type(v, c, r); }
};

template<typename iterator_type>
struct base_vector_iterator
{
    size_t index;
    explicit base_vector_iterator(size_t i) : index(i) {}
    iterator_type &operator++() { ++index; return *static_cast<iterator_type*>(this); }
private:
    base_vector_iterator &operator=(const base_vector_iterator&);
};

template<typename t>
struct iterator
{
    struct type : public base_vector_iterator<type>
    {
        t &ref;
        type(t &v, size_t i) : base_vector_iterator(i), ref(v) { ASSERT(index<ref.size()); }
        typename t::reference operator*() { ASSERT(index<ref.size()); return ref(index); }
    };
    static type begin(t &v, size_t i) { return type(v, i); }
};

template<typename t>
struct const_iterator
{
    struct type : public base_vector_iterator<type>
    {
        const t &ref;
        type(const t &v, size_t i) : base_vector_iterator(i), ref(v) { ASSERT(index<ref.size()); }
        typename t::const_reference operator*() const { ASSERT(index<ref.size()); return ref(index); }
    };
    static type begin(const t &v, size_t i) { return type(v, i); }
};


template<typename value_t, typename expr_t, typename tag> inline
typename iterator<expr_t>::type begin(vector_expression<value_t, expr_t, tag> &v, size_t i = 0)
{ return iterator<expr_t>::begin(get_vector(v), i); }

template<typename value_t, typename expr_t, typename tag> inline
typename const_iterator<expr_t>::type begin(const vector_expression<value_t, expr_t, tag> &v, size_t i = 0)
{ return const_iterator<expr_t>::begin(get_vector(v), i); }


template<typename value_t, typename expr_t, typename tag_t>
class column_wrapper :
    public matrix_expression<value_t,
                             column_wrapper<value_t, expr_t, tag_t>,
                             tag::expression<util::tag_complexity<tag_t>::value>
                             >
{
public:
    typedef value_t value_type; 
    typedef const value_type const_reference;   

    column_wrapper()
    {}

    const_reference operator()(size_t row, size_t col) const
    {
        ASSERT(col==0);
        return (*static_cast<const expr_t*>(this))(row);
    }

    size_t rows() const { return static_cast<const expr_t*>(this)->size(); }

    size_t cols() const { return 1; }
private:
    column_wrapper(const column_wrapper&);
    column_wrapper &operator=(const column_wrapper&);
};


template<typename value_t, typename expr_t, typename tag>
class vector_expression<value_t, expr_t, tag, orientation::column> : public column_wrapper<value_t, expr_t, tag>
{
private:    //hiding members of 'column_wrapper'
    using column_wrapper<value_t, expr_t, tag>::value_type;
    using column_wrapper<value_t, expr_t, tag>::const_reference;
    using column_wrapper<value_t, expr_t, tag>::operator();
    using column_wrapper<value_t, expr_t, tag>::rows;
    using column_wrapper<value_t, expr_t, tag>::cols;

    friend orientation::column orientation_type(expr_t*) { return orientation::column(); }
};


template<typename value_t, typename expr_t, typename tag_t>
class row_wrapper :
    public matrix_expression<value_t,
                             row_wrapper<value_t, expr_t, tag_t>,
                             tag::expression<util::tag_complexity<tag_t>::value>
                             >
{
public:
    typedef value_t value_type; 
    typedef const value_type const_reference;   

    row_wrapper()
    {}

    const_reference operator()(size_t row, size_t col) const
    {
        ASSERT(row==0);
        return (*static_cast<const expr_t*>(this))(col);
    }

    size_t rows() const { return 1; }

    size_t cols() const { return static_cast<const expr_t*>(this)->size(); }
private:
    row_wrapper(const row_wrapper&);
    row_wrapper &operator=(const row_wrapper&);
};


template<typename value_t, typename expr_t, typename tag>
class vector_expression<value_t, expr_t, tag, orientation::row> : public row_wrapper<value_t, expr_t, tag>
{
private:    //hiding members of 'row_wrapper'
    using row_wrapper<value_t, expr_t, tag>::value_type;
    using row_wrapper<value_t, expr_t, tag>::const_reference;
    using row_wrapper<value_t, expr_t, tag>::operator();
    using row_wrapper<value_t, expr_t, tag>::rows;
    using row_wrapper<value_t, expr_t, tag>::cols;

    friend orientation::row orientation_type(expr_t*) { return orientation::row(); }
};

template<typename value_t, typename expr_t, typename tag_t, typename orientation_t>
struct expr_traits<vector_expression<value_t, expr_t, tag_t, orientation_t> >
{
    typedef value_t value_type; 
    typedef expr_t type;    
    typedef tag_t tag;  
    typedef orientation_t orientation;  
    enum {
        complexity = util::tag_complexity<tag>::value   
    };
};




template<typename value_t, typename expr_t, typename tag> inline
expr_t &get_vector(vector_expression<value_t, expr_t, tag> &a)
{ return *static_cast<expr_t*>(&a); }


template<typename value_t, typename expr_t, typename tag> inline
const expr_t &get_vector(const vector_expression<value_t, expr_t, tag> &a)
{ return *static_cast<const expr_t*>(&a); }

namespace util {
template<typename type, size_t static_size>
class sequence_base
{
public:
    sequence_base() : _p(0) {}
    ~sequence_base() { if (_p && _p!=_data) delete[] _p; }
    void realloc(size_t size)
    {
        type *p = size<=static_size ? _data : new type[size];   //if (size>static_size) allocates memory
        if (_p && _p!=_data) delete[] _p;                       //note that static_size>0
        _p = size!=0 ? p : 0;   //if size==0 frees memory
        return;
    }
protected:
    type *_p;
    type _data[static_size];
};
template<typename type>
class sequence_base<type, 0>
{
public:
    sequence_base() : _p(0) {}
    ~sequence_base() { if (_p) delete[] _p; }
    void realloc(size_t size)
    {
        type *p = size!=0 ? new type[size] : 0; //if size==0 frees memory
        if (_p) delete[] _p;
        _p = p;
        return;
    }
protected:
    type *_p;
};
template<typename type, size_t static_size = 0> class sequence;
template<typename type, size_t static_size> void swap(sequence<type, static_size>&, size_t, sequence<type, static_size>&, size_t);
template<typename type> void swap(sequence<type>&, sequence<type>&);
template<typename type, size_t static_size>
class sequence : public sequence_base<type, static_size>
{
public:
    typedef type value_type;    
    typedef type &reference;    
    typedef const type &const_reference;    

    sequence() {}
    sequence(size_t size) { this->realloc(size); }
    sequence(size_t size, const value_type &value) { this->realloc(size, value); }
    value_type &operator[](size_t i) { ASSERT(!this->is_null()); return _p[i]; }
    const value_type &operator[](size_t i) const { ASSERT(!this->is_null()); return _p[i]; }
    template<size_t s>
    sequence &copy(const sequence<type, s> &a, size_t size)
    {
        ASSERT(!this->is_null() && !a.is_null());
        for (size_t i = 0; i<size; ++i)
            _p[i] = a._p[i];
        return *this;
    }
    void fill(const value_type &value, size_t count)
    {
        ASSERT(!this->is_null());
        for (size_t i = 0; i<count; ++i)
            _p[i] = value;
        return;
    }
    using sequence_base<type, static_size>::realloc;
    void realloc(size_t size, const value_type &value)
    {
        this->realloc(size);
        fill(value, size);
        return;
    }
    bool is_null() const { return _p==0; }

private:
    //value_type *_p;
    //value_type _data[static_size];

    sequence(const sequence&);
    sequence &operator=(const sequence&);

    friend void swap<type, static_size>(sequence<type, static_size>&, size_t, sequence<type, static_size>&, size_t);
    friend void swap<type>(sequence<type>&, sequence<type>&);

};
template<typename type, size_t static_size>
void swap(sequence<type, static_size> &a, size_t a_size, sequence<type, static_size> &b, size_t b_size)
{
    using std::swap;
    const bool
        a_is_static = a_size<=static_size,
        b_is_static = b_size<=static_size;
    if (!a_is_static)
    {
        if (!b_is_static)
            swap(a._p, b._p);
        else
        {
            for (size_t i = 0; i<b_size; ++i)
                a._data[i] = b._p[i];
            b._p = a._p;
            a._p = a._data;
        }
        return;
    }
    else if (!b_is_static)
    {
        for (size_t i = 0; i<a_size; ++i)
            b._data[i] = a._p[i];
        a._p = b._p;
        b._p = b._data;
        return;
    }
    size_t min;
    if (a_size<b_size)
    {
        min = a_size;
        for (size_t i = min; i<b_size; ++i)
            a._p[i] = b._p[i];
    }
    else//a_size>=b_size
    {
        min = b_size;
        for (size_t i = min; i<a_size; ++i)
            b._p[i] = a._p[i];
    }
    for (size_t i = 0; i<min; ++i)
        swap(a._p[i], b._p[i]);
    return;
}
template<typename type> inline
void swap(sequence<type> &a, sequence<type> &b)
{
    using std::swap;
    swap(a._p, b._p);
    return;
}
}//namespace util


template<typename type>
class temporary : public type
{
public:
    temporary() : type() {}

    temporary(temporary &a) : type()
    {
        using std::swap;
        swap(*static_cast<type*>(this), *static_cast<type*>(&a));
    }

    template<typename other>
    temporary(const other &a) : type(a) {}

    temporary &operator=(temporary &a)
    {
        using std::swap;
        swap(*static_cast<type*>(this), *static_cast<type*>(&a));
        return *this;
    }

    template<typename other>
    temporary &operator=(const other &a) { type::operator=(a); return *this; }
};


template<typename type>
struct rvalue : public type
{
    rvalue();   //these member declarations are needed for compiler not to complain about the impossibility
    ~rvalue();  //to generate corresponding members (these two shall not be defined!!!)
};




template<typename type>
rvalue<type> &to_rvalue_ref(type &a)
{ return *static_cast<rvalue<type>*>(&a); }

//forward declaration
template<typename expression_type, typename tag, typename orientation> class vector_ref;


template<typename expression_type>
class vector_ref<expression_type, tag::reference, orientation::row> :
    public vector_expression<typename expr_traits<expression_type>::value_type,
                             vector_ref<expression_type, tag::reference, orientation::row>,
                             tag::reference,
                             orientation::row
                             >
{
    typedef typename expr_traits<expression_type>::type expr_t;
public:
    typedef typename expr_t::value_type value_type; 
    typedef typename expr_t::reference reference;   
    typedef typename expr_t::const_reference const_reference;   

    vector_ref(vector_ref &a) : _ref(a._ref), _row(a._row)
    {}
    vector_ref(rvalue<vector_ref> &a) : _ref(a._ref), _row(a._row)
    {}
    vector_ref(expr_t &a, size_t row) : _ref(a), _row(row)
    {}
    reference operator()(size_t i)
    { return _ref(_row, i); }
    const_reference operator()(size_t i) const
    { return _ref(_row, i); }
    reference operator[](size_t i)
    { return _ref(_row, i); }
    const_reference operator[](size_t i) const
    { return _ref(_row, i); }
    size_t size() const
    { return _ref.cols(); }
    //Conversion to "reference to rvalue".
    operator rvalue<vector_ref>&()
    { return to_rvalue_ref(*this); }

private:
    expr_t &_ref;
    const size_t _row;

    vector_ref &operator=(const vector_ref&);   //do not implement

    friend expr_t &get_vector_ref(vector_ref &a) { return a._ref; }
    friend const expr_t &get_vector_ref(const vector_ref &a) { return a._ref; }
    friend size_t get_vector_row(const vector_ref &a) { return a._row; }
};

namespace util {
template<size_t a, size_t b, size_t c = 0>
struct select
{
    template<bool cond, size_t true_clause, size_t false_clause>
    struct _if { enum { value = true_clause }; };
    template<size_t true_clause, size_t false_clause>
    struct _if<false, true_clause, false_clause> { enum { value = false_clause }; };
    enum {
        max_of_2 = _if<(a>b), a, b>::value, 
        max = _if<(max_of_2>c), max_of_2, c>::value 
    };
};
}//namespace util


template<typename expression_type, size_t complexity>
class vector_ref<expression_type, tag::expression<complexity>, orientation::row> :
    public vector_expression<typename expr_traits<expression_type>::value_type,
                             vector_ref<expression_type, tag::expression<complexity>, orientation::row>,
                             tag::expression<util::select<complexity, 1>::max - 1>,
                             orientation::row
                             >
{
    typedef typename expr_traits<expression_type>::type expr_t;
public:
    typedef typename expr_t::value_type value_type; 
    typedef typename expr_t::const_reference const_reference;   

    vector_ref(const vector_ref &a) : _ref(a._ref), _row(a._row)
    {}
    vector_ref(const expr_t &a, size_t row) : _ref(a), _row(row)
    {}
    const_reference operator()(size_t i) const
    { return _ref(_row, i); }
    const_reference operator[](size_t i) const
    { return _ref(_row, i); }
    size_t size() const
    { return _ref.cols(); }

private:
    const expr_t &_ref;
    const size_t _row;

    vector_ref &operator=(const vector_ref&);   //do not implement

    friend expr_t &get_vector_ref(vector_ref &a) { return a._ref; }
    friend const expr_t &get_vector_ref(const vector_ref &a) { return a._ref; }
    friend size_t get_vector_row(const vector_ref &a) { return a._row; }
};



template<typename expression_type, typename tag> inline
typename util::non_zero_in_row<typename expr_traits<expression_type>::shape>::first_type
first_non_zero(const vector_ref<expression_type, tag, orientation::row> &v)
{ return first_non_zero_in_row(get_vector_ref(v), get_vector_row(v)); }
template<typename expression_type, typename tag> inline
typename util::non_zero_in_row<typename expr_traits<expression_type>::shape>::last_type
last_non_zero(const vector_ref<expression_type, tag, orientation::row> &v)
{ return last_non_zero_in_row(get_vector_ref(v), get_vector_row(v)); }


template<typename expression_type, typename tag>
struct iterator<vector_ref<expression_type, tag, orientation::row> >
{
    typedef typename row_iterator<typename expr_traits<expression_type>::type>::type type;
    static type begin(vector_ref<expression_type, tag, orientation::row> &v, size_t i = 0)
    { return row_iterator<typename expr_traits<expression_type>::type>::begin(get_vector_ref(v), get_vector_row(v), i); }
};

template<typename expression_type, typename tag>
struct const_iterator<vector_ref<expression_type, tag, orientation::row> >
{
    typedef typename row_const_iterator<typename expr_traits<expression_type>::type>::type type;
    static type begin(const vector_ref<expression_type, tag, orientation::row> &v, size_t i = 0)
    { return row_const_iterator<typename expr_traits<expression_type>::type>::begin(get_vector_ref(v), get_vector_row(v), i); }
};


struct access_error : public std::logic_error
{
    access_error() : std::logic_error("access violation (inaccessible element)") {}
    virtual ~access_error() throw() {}
};

namespace util {
struct details_base
{
    static size_t alloc_size(size_t size) { return (size*(size + 1))/2; }
};
//forward declaration
template<typename shape> struct details;
template<>
struct details<lower_triangular> : public details_base
{
    static size_t index(size_t row, size_t column, size_t/*size*/ = 0)
    { return details_base::alloc_size(row) + column; }
    template<typename storage>
    static typename storage::reference element(storage &p, size_t row, size_t column, size_t/*size*/)
    {
        if (row<column)
            throw access_error();
        return p[details::index(row, column)];
    }
    template<typename storage>
    static typename storage::const_reference element(const storage &p, size_t row, size_t column, size_t/*size*/)
    { return row>=column ? p[details::index(row, column)] : traits<typename storage::value_type>::zero; }
    template<typename storage, typename value_t, typename expr_t, typename tag>
    static void copy_values(storage &p, const shaped_expression<value_t, expr_t, tag, lower_triangular> &rhs, size_t size)
    {
        const expr_t &a = get_matrix(rhs);
        for (size_t i = 0, k = 0; i<size; ++i)
            for (size_t j = 0; j<=i; ++j, ++k)
                p[k] = a(i, j);
        return;
    }
};
template<>
struct details<symmetric> : public details_base
{
    static size_t index(size_t row, size_t column, size_t/*size*/ = 0)
    {
        typedef details<lower_triangular> det_lower_tr;
        return row>=column ? det_lower_tr::index(row, column) : det_lower_tr::index(column, row);
    }
    template<typename storage>
    static typename storage::reference element(storage &p, size_t row, size_t column, size_t/*size*/)
    { return p[details::index(row, column)]; }
    template<typename storage>
    static typename storage::const_reference element(const storage &p, size_t row, size_t column, size_t/*size*/)
    { return p[details::index(row, column)]; }
    template<typename storage, typename value_t, typename expr_t, typename tag>
    static void copy_values(storage &p, const shaped_expression<value_t, expr_t, tag, symmetric> &rhs, size_t size)
    {
        const expr_t &a = get_matrix(rhs);
        for (size_t i = 0, k = 0; i<size; ++i)
            for (size_t j = 0; j<=i; ++j, ++k)
                p[k] = a(i, j);
        return;
    }
};
template<>
struct details<upper_triangular> : public details_base
{
    static size_t index(size_t row, size_t column, size_t size)
    { return (row*(2*size - row - 1))/2 + column; } //returns n'th element in triangular matrix stored row-by-row
    template<typename storage>
    static typename storage::reference element(storage &p, size_t row, size_t column, size_t size)
    {
        if (row>column)
            throw access_error();
        return p[details::index(row, column, size)];
    }
    template<typename storage>
    static typename storage::const_reference element(const storage &p, size_t row, size_t column, size_t size)
    { return row<=column ? p[details::index(row, column, size)] : traits<typename storage::value_type>::zero; }
    template<typename storage, typename value_t, typename expr_t, typename tag>
    static void copy_values(storage &p, const shaped_expression<value_t, expr_t, tag, upper_triangular> &rhs, size_t size)
    {
        const expr_t &a = get_matrix(rhs);
        for (size_t i = 0, k = 0; i<size; ++i)
            for (size_t j = i; j<size; ++j, ++k)
                p[k] = a(i, j);
        return;
    }
};
}//namespace util

//forward declaration
template<typename value_t, typename shape = rectangular> class matrix;


template<typename value_t, typename shape>
class matrix : public shaped_expression<value_t, matrix<value_t, shape>, tag::storage, shape>
{
    typedef shaped_expression<value_t, matrix<value_t, shape>, tag::storage, shape> expression_type;
public:
    typedef value_t value_type; 
    typedef typename util::sequence<value_type>::reference reference;   
    typedef typename util::sequence<value_type>::const_reference const_reference;   
    typedef vector_ref<expression_type, tag::reference, orientation::row> row_ref;
    typedef vector_ref<expression_type, tag::expression<0>, orientation::row> const_row_ref;

    matrix() : _p(), _size(0)
    {}
    matrix(const matrix &a) : _p(), _size(0)
    { this->operator=(a); }
    matrix(temporary<matrix> &a) : _p(), _size(0)
    { this->swap(a); }
    template<typename expr_t, typename tag>
    matrix(const shaped_expression<value_type, expr_t, tag, shape> &a) : _p(), _size(0)
    { this->assign(a); }
    explicit matrix(size_t size) : _p(), _size(0)
    { this->resize(size); }
    matrix(size_t size, const_reference value) : _p(), _size(0)
    { this->resize(size, value); }
    ~matrix()
    {}
    matrix &operator=(const matrix &a)
    {
        this->resize(a._size);
        _p.copy(a._p, util::details<shape>::alloc_size(_size));
        return *this;
    }
    template<typename expr_t, typename tag>
    matrix &operator=(const shaped_expression<value_type, expr_t, tag, shape> &a)
    {
        matrix(a).swap(*this);
        return *this;
    }
    matrix &operator=(temporary<matrix> &a)
    { this->swap(a); return *this; }

    reference operator()(size_t row, size_t column)
    {
        ASSERT(row<this->rows() && column<this->cols());
        return util::details<shape>::element(_p, row, column, _size);
    }
    const_reference operator()(size_t row, size_t column) const
    {
        ASSERT(row<this->rows() && column<this->cols());
        return util::details<shape>::element(_p, row, column, _size);
    }
    row_ref operator()(size_t row)
    { ASSERT(row<this->rows()); return row_ref(*this, row); }
    const_row_ref operator()(size_t row) const
    { ASSERT(row<this->rows()); return const_row_ref(*this, row); }
    row_ref operator[](size_t row)
    { ASSERT(row<this->rows()); return row_ref(*this, row); }
    const_row_ref operator[](size_t row) const
    { ASSERT(row<this->rows()); return const row_ref(*this, row); }

    template<typename expr_t, typename tag>
    void assign(const shaped_expression<value_type, expr_t, tag, shape> &rhs)
    {
        const expr_t &a = get_matrix(rhs);
        ASSERT(a.rows()==a.cols()); //equal by definition
        this->resize(a.rows());
        util::details<shape>::copy_values(_p, rhs, _size);
        return;
    }

    void resize(size_t size)
    {
        if (size!=_size)
        {
            _p.realloc(util::details<shape>::alloc_size(size));
            _size = size;
        }
        return;
    }

    void resize(size_t size, const_reference value)
    {
        this->resize(size);
        _p.fill(value, util::details<shape>::alloc_size(_size));
        return;
    }

    void swap(matrix &a)
    {
        using std::swap;
        swap(_p, a._p);
        swap(_size, a._size);
        return;
    }
    size_t rows() const
    { return _size; }
    size_t cols() const
    { return _size; }

private:
    util::sequence<value_type> _p;
    size_t _size;
};

namespace util {
struct base_symmetric_matrix_iterator
{
    const size_t row;
    size_t col;
    base_symmetric_matrix_iterator(size_t r, size_t c) : row(r), col(c) {}
    template<typename ptr_t>
    void advance(ptr_t &ptr)
    {
        ptr += size_t(row<=col)*col + 1; ++col;
        return;
    }
private:
    base_symmetric_matrix_iterator &operator=(const base_symmetric_matrix_iterator&);
};
template<typename value_t, typename iterator_type>
struct symmetric_matrix_iterator : public base_symmetric_matrix_iterator
{
    typename matrix<value_t, symmetric>::value_type *ptr;
    symmetric_matrix_iterator(matrix<value_t, symmetric> &m, size_t r, size_t c) :
        base_symmetric_matrix_iterator(r, c), ptr(&m(r, c))
    {}
    typename matrix<value_t, symmetric>::reference operator*() { return *ptr; }
    iterator_type &operator++() { this->advance(ptr); return *static_cast<iterator_type*>(this); }
};
template<typename value_t, typename iterator_type>
struct symmetric_matrix_const_iterator : public base_symmetric_matrix_iterator
{
    const typename matrix<value_t, symmetric>::value_type *ptr;
    symmetric_matrix_const_iterator(const matrix<value_t, symmetric> &m, size_t r, size_t c) :
        base_symmetric_matrix_iterator(r, c), ptr(&m(r, c))
    {}
    typename matrix<value_t, symmetric>::const_reference operator*() const { return *ptr; }
    iterator_type &operator++() { this->advance(ptr); return *static_cast<iterator_type*>(this); }
};
}//namespace util

template<typename value_t>
struct row_iterator<matrix<value_t, symmetric> >
{
    struct type : public util::symmetric_matrix_iterator<value_t, type>
    {
        type(matrix<value_t> &m, size_t r, size_t c) : util::symmetric_matrix_iterator(m, r, c) {}
    };
    static type begin(matrix<value_t, symmetric> &m, size_t r, size_t c) { return type(m, r, c); }
};
template<typename value_t>
struct row_const_iterator<matrix<value_t, symmetric> >
{
    struct type : public util::symmetric_matrix_const_iterator<value_t, type>
    {
        type(const matrix<value_t, symmetric> &m, size_t r, size_t c) : symmetric_matrix_const_iterator(m, r, c) {}
    };
    static type begin(const matrix<value_t, symmetric> &m, size_t r, size_t c) { return type(m, r, c); }
};

template<typename value_t>
struct column_iterator<matrix<value_t, symmetric> >
{
    struct type : public util::symmetric_matrix_iterator<value_t, type>
    {
        type(matrix<value_t, symmetric> &m, size_t c, size_t r) : util::symmetric_matrix_iterator(m, c, r) {}
    };
    static type begin(matrix<value_t, symmetric> &m, size_t c, size_t r) { return type(m, c, r); }
};
template<typename value_t>
struct column_const_iterator<matrix<value_t, symmetric> >
{
    struct type : public util::symmetric_matrix_const_iterator<value_t, type>
    {
        type(const matrix<value_t, symmetric> &m, size_t c, size_t r) : symmetric_matrix_const_iterator(m, c, r) {}
    };
    static type begin(const matrix<value_t, symmetric> &m, size_t c, size_t r) { return type(m, c, r); }
};

//forward declaration
template<typename expression1, typename expression2> class product;


template<typename value_t>
class matrix<value_t, rectangular> :
    public matrix_expression<value_t, matrix<value_t, rectangular>, tag::storage>
{
    typedef matrix_expression<value_t, matrix<value_t, rectangular>, tag::storage> expression_type;
public:
    typedef value_t value_type; 
    typedef typename util::sequence<value_type>::reference reference;   
    typedef typename util::sequence<value_type>::const_reference const_reference;   
    typedef vector_ref<expression_type, tag::reference, orientation::row> row_ref;
    typedef vector_ref<expression_type, tag::expression<0>, orientation::row> const_row_ref;


    matrix() : _p(), _rows(0), _cols(0)
    {}
    matrix(const matrix &a) : _p(), _rows(0), _cols(0)
    { operator=(a); }
    matrix(temporary<matrix> &a) : _p(), _rows(0), _cols(0)
    { this->swap(a); }
    template<typename expr_t, typename tag>
    matrix(const matrix_expression<value_type, expr_t, tag> &a) : _p(), _rows(0), _cols(0)
    { this->assign(a); }
    matrix(size_t nrows, size_t ncols) : _p(), _rows(0), _cols(0)
    { this->resize(nrows, ncols); }
    matrix(size_t nrows, size_t ncols, const_reference value) : _p(), _rows(0), _cols(0)
    { this->resize(nrows, ncols, value); }
    ~matrix()
    {}
    matrix &operator=(const matrix &a)
    {
        this->resize(a._rows, a._cols);
        _p.copy(a._p, _rows*_cols);
        return *this;
    }
    template<typename expr_t, typename tag>
    matrix &operator=(const matrix_expression<value_type, expr_t, tag> &rhs)
    {
        matrix(rhs).swap(*this);
        return *this;
    }
    matrix &operator=(temporary<matrix> &a)
    { this->swap(a); return *this; }
    reference operator()(size_t row, size_t column)
    { ASSERT(row<this->rows() && column<this->cols()); return _p[row*_cols + column]; }
    const_reference operator()(size_t row, size_t column) const
    { ASSERT(row<this->rows() && column<this->cols()); return _p[row*_cols + column]; }
    row_ref operator()(size_t row)
    { ASSERT(row<this->rows()); return row_ref(*this, row); }
    const_row_ref operator()(size_t row) const
    { ASSERT(row<this->rows()); return const_row_ref(*this, row); }
    row_ref operator[](size_t row)
    { ASSERT(row<this->rows()); return row_ref(*this, row); }
    const_row_ref operator[](size_t row) const
    { ASSERT(row<this->rows()); return const_row_ref(*this, row); }

    template<typename expr_t, typename tag>
    void assign(const matrix_expression<value_type, expr_t, tag> &rhs)
    {
        const expr_t &a = get_matrix(rhs);
        this->resize(a.rows(), a.cols());
        for (size_t i = 0, k = 0; i<_rows; ++i)
            for (size_t j = 0; j<_cols; ++j, ++k)
                _p[k] = a(i, j);
        return;
    }
    template<typename expression1, typename expression2, typename tag>
    void assign(const matrix_expression<value_type, product<expression1, expression2>, tag> &rhs);

    void resize(size_t nrows, size_t ncols)
    {
        const size_t size = nrows*ncols;
        if (size==0)
            nrows = ncols = 0;
        if (size!=_rows*_cols)
        {
            _p.realloc(size);
            _rows = nrows;
            _cols = ncols;
        }
        return;
    }

    void resize(size_t nrows, size_t ncols, const_reference value)
    {
        this->resize(nrows, ncols);
        _p.fill(value, _rows*_cols);
        return;
    }

    void swap(matrix &a)
    {
        using std::swap;
        swap(_p, a._p);
        swap(_rows, a._rows);
        swap(_cols, a._cols);
        return;
    }
    size_t rows() const
    { return _rows; }
    size_t cols() const
    { return _cols; }

private:
    util::sequence<value_type> _p;
    size_t _rows, _cols;
};

template<typename value_t>
struct row_iterator<matrix<value_t> >
{
    struct type
    {
        typename matrix<value_t>::value_type *ptr;
        type(matrix<value_t> &m, size_t r, size_t c) : ptr(&m(r, c)) {}
        type &operator++() { ++ptr; return *this; }
        typename matrix<value_t>::reference operator*() { return *ptr; }
    private:
        type &operator=(const type&);
    };
    static type begin(matrix<value_t> &m, size_t r, size_t c) { return type(m, r, c); }
};
template<typename value_t>
struct row_const_iterator<matrix<value_t> >
{
    struct type
    {
        const typename matrix<value_t>::value_type *ptr;
        type(const matrix<value_t> &m, size_t r, size_t c) : ptr(&m(r, c)) {}
        type &operator++() { ++ptr; return *this; }
        typename matrix<value_t>::const_reference operator*() const { return *ptr; }
    private:
        type &operator=(const type&);
    };
    static type begin(const matrix<value_t> &m, size_t r, size_t c) { return type(m, r, c); }
};
template<typename value_t>
struct column_iterator<matrix<value_t> >
{
    struct type
    {
        typename matrix<value_t>::value_type *ptr;
        const size_t cols;
        type(matrix<value_t> &m, size_t c, size_t r) : ptr(&m(r, c)), cols(m.cols()) {}
        type &operator++() { ptr += cols; return *this; }
        typename matrix<value_t>::reference operator*() { return *ptr; }
    private:
        type &operator=(const type&);
    };
    static type begin(matrix<value_t> &m, size_t c, size_t r) { return type(m, c, r); }
};
template<typename value_t>
struct column_const_iterator<matrix<value_t> >
{
    struct type
    {
        const typename matrix<value_t>::value_type *ptr;
        const size_t cols;
        type(const matrix<value_t> &m, size_t c, size_t r) : ptr(&m(r, c)), cols(m.cols()) {}
        type &operator++() { ptr += cols; return *this; }
        typename matrix<value_t>::const_reference operator*() const { return *ptr; }
    private:
        type &operator=(const type&);
    };
    static type begin(const matrix<value_t> &m, size_t c, size_t r) { return type(m, c, r); }
};


template<typename value_t>
class vector : public vector_expression<value_t, vector<value_t>, tag::storage, orientation::column>
{
public:
    typedef value_t value_type; 
    typedef typename util::sequence<value_type>::reference reference;   
    typedef typename util::sequence<value_type>::const_reference const_reference;   


    vector() : _p(), _size(0)
    {}
    vector(const vector &a) : _p(), _size(0)
    { operator=(a); }
    vector(temporary<vector> &a) : _p(), _size(0)
    { this->swap(a); }
    template<typename expr_t, typename tag>
    vector(const vector_expression<value_type, expr_t, tag> &a) : _p(), _size(0)
    { this->assign(a); }
    explicit vector(size_t size) : _p(), _size(0)
    { this->resize(size); }
    vector(size_t size, const_reference value) : _p(), _size(0)
    { this->resize(size, value); }
    ~vector()
    {}
    vector &operator=(const vector &a)
    {
        this->resize(a._size);
        _p.copy(a._p, _size);
        return *this;
    }
    template<typename expr_t, typename tag>
    vector &operator=(const vector_expression<value_type, expr_t, tag> &rhs)
    {
        vector<value_type>(rhs).swap(*this);
        return *this;
    }
    vector &operator=(temporary<vector> &a)
    { this->swap(a); return *this; }
    reference operator()(size_t i)
    { ASSERT(i<this->size()); return _p[i]; }
    const_reference operator()(size_t i) const
    { ASSERT(i<this->size()); return _p[i]; }
    reference operator[](size_t i)
    { ASSERT(i<this->size()); return _p[i]; }
    const_reference operator[](size_t i) const
    { ASSERT(i<this->size()); return _p[i]; }

    template<typename expr_t, typename tag>
    void assign(const vector_expression<value_type, expr_t, tag> &rhs)
    {
        const expr_t &a = get_vector(rhs);
        this->resize(a.size());
        for (size_t i = 0; i<_size; ++i)
            _p[i] = a(i);
        return;
    }

    void resize(size_t size)
    {
        if (size!=_size)
        {
            _p.realloc(size);
            _size = size;
        }
        return;
    }

    void resize(size_t size, const_reference value)
    {
        this->resize(size);
        _p.fill(value, _size);
        return;
    }

    void swap(vector &a)
    {
        using std::swap;
        swap(_p, a._p);
        swap(_size, a._size);
        return;
    }
    size_t size() const
    { return _size; }

private:
    util::sequence<value_type> _p;
    size_t _size;
};




template<typename value_t, typename expr_t> inline
void swap(matrix_expression<value_t, expr_t, tag::storage> &a,
          matrix_expression<value_t, expr_t, tag::storage> &b)
{ get_matrix(a).swap(get_matrix(b)); return; }


template<typename value_t, typename expr_t, typename shape> inline
void swap(shaped_expression<value_t, expr_t, tag::storage, shape> &a,
          shaped_expression<value_t, expr_t, tag::storage, shape> &b)
{ get_matrix(a).swap(get_matrix(b)); return; }


struct reshape_error : public std::logic_error
{
    reshape_error() : std::logic_error("matrix must be square (of size strictly n by n)") {}
    virtual ~reshape_error() throw() {}
};


struct size_error : public std::invalid_argument
{
    explicit size_error(const std::string &msg) : std::invalid_argument(msg) {}
    virtual ~size_error() throw() {}
};


template<typename value_t, typename expr_t, typename tag, typename shape>
struct base_matrix_expression :
    public shaped_expression<value_t, expr_t, tag, shape>
{};

template<typename value_t, typename expr_t, typename tag>
struct base_matrix_expression<value_t, expr_t, tag, rectangular> :
    public matrix_expression<value_t, expr_t, tag>
{};

namespace util {
template<typename side> struct choose;
template<>
struct choose<upper_triangular>
{
    template<typename expr_t>
    static typename expr_t::reference element(expr_t &a, size_t row, size_t col)
    { return a(row, col); }
    template<typename expr_t>
    static typename expr_t::const_reference element(const expr_t &a, size_t row, size_t col)
    { return a(row, col); }
};
template<>
struct choose<lower_triangular>
{
    template<typename expr_t>
    static typename expr_t::reference element(expr_t &a, size_t row, size_t col)
    { return a(col, row); }
    template<typename expr_t>
    static typename expr_t::const_reference element(const expr_t &a, size_t row, size_t col)
    { return a(col, row); }
};
template<typename shape, typename side = upper_triangular>
struct translate_access
{
    template<typename expr_t>
    static typename expr_t::reference get(expr_t &a, size_t row, size_t col)
    { return choose<side>::element(a, row, col); }
    template<typename expr_t>
    static typename expr_t::const_reference get(const expr_t &a, size_t row, size_t col)
    { return choose<side>::element(a, row, col); }
};
//struct minmax
//{
//  const size_t x, y, mask;
//  minmax(size_t a, size_t b) : x(a), y(b), mask((x ^ y) & -(x < y))
//  { ASSERT(std::min(a, b)==min() && std::max(a, b)==max()); }
//  size_t min() const { return y^mask; }
//  size_t max() const { return x^mask; }
//};
struct minmax
{
    size_t sup, inf;
    minmax(size_t a, size_t b)
    {
        if (a<b)
        {
            sup = a;
            inf = b;
        }
        else
        {
            sup = b;
            inf = a;
        }
        ASSERT(std::min(a, b)==min() && std::max(a, b)==max());
    }
    size_t min() const { return sup; }
    size_t max() const { return inf; }
};
template<typename side>
struct translate_access<symmetric, side>
{
    template<typename expr_t>
    static typename expr_t::reference get(expr_t &a, size_t row, size_t col)
    {
        //const size_t
        //  r = std::min(row, col),
        //  c = std::max(row, col);
        //return choose<side>::element(a, r, c);
        const minmax i(row, col);
        return choose<side>::element(a, i.min(), i.max());
    }
    template<typename expr_t>
    static typename expr_t::const_reference get(const expr_t &a, size_t row, size_t col)
    {
        //const size_t
        //  r = std::min(row, col),
        //  c = std::max(row, col);
        //return choose<side>::element(a, r, c);
        const minmax i(row, col);
        return choose<side>::element(a, i.min(), i.max());
    }
};
template<typename side>
struct translate_access<upper_triangular, side>
{

    template<typename expr_t>
    static typename expr_t::reference get(expr_t &a, size_t row, size_t col)
    {
        if (row<=col)
            return choose<side>::element(a, row, col);
        throw access_error();
    }
    template<typename expr_t>
    static typename expr_t::const_reference get(const expr_t &a, size_t row, size_t col)
    { return row<=col ? choose<side>::element(a, row, col) : traits<typename expr_t::value_type>::zero; }
};
template<typename side>
struct translate_access<lower_triangular, side>
{

    template<typename expr_t>
    static typename expr_t::reference get(expr_t &a, size_t row, size_t col)
    {
        if (row>=col)
            return choose<side>::element(a, col, row);
        throw access_error();
    }
    template<typename expr_t>
    static typename expr_t::const_reference get(const expr_t &a, size_t row, size_t col)
    { return row>=col ? choose<side>::element(a, col, row) : traits<typename expr_t::value_type>::zero; }
};
template<typename expr_t>
struct const_ref
{
    const expr_t &ref;  

    explicit const_ref(const expr_t &a) : ref(a)
    {}
    operator const expr_t&() { return ref; }

private:
    const_ref &operator=(const const_ref&); //do not define
};
template<typename expr_t> inline
const_ref<expr_t> make_const_ref(const expr_t &a)
{ return const_ref<expr_t>(a); }
}//namespace util


template<typename target_shape>
struct reshape_constraint
{

    template<typename value_t, typename expr_t, typename tag>
    static bool test(const matrix_expression<value_t, expr_t, tag> &a)
    {
        const expr_t &m = get_matrix(a);
        if (m.rows()!=m.cols()) //matrix must be square!
            throw reshape_error();
        return true;
    }
};

template<>
struct reshape_constraint<rectangular>
{

    template<typename value_t, typename expr_t, typename tag>
    static bool test(const matrix_expression<value_t, expr_t, tag>&)
    { return true; }
};

template<typename expression_type, typename target_shape, typename side = upper_triangular> class reshaped;

template<typename value_t, typename expr_t, typename tag_t, typename target_shape>
class reshaped<matrix_expression<value_t, expr_t, tag_t>, target_shape, upper_triangular> :
    public base_matrix_expression<
            typename expr_t::value_type,
            reshaped<matrix_expression<value_t, expr_t, tag_t>, target_shape, upper_triangular>,
            tag::expression<util::tag_complexity<tag_t>::value>,
            target_shape
                                  >
{
public:
    typedef typename expr_t::value_type value_type; 
    //typedef typename expr_t::reference reference; ///<Type of reference to element.
    typedef typename expr_t::const_reference const_reference;   

    reshaped(const util::const_ref<expr_t> &a) : _ref(a.ref)
    {   //assertion is needed if user tries to construct reshaped object directly
        ASSERT(reshape_constraint<target_shape>::test(a.ref));
    }
    const_reference operator()(size_t row, size_t col) const
    { return util::translate_access<target_shape, upper_triangular>::get(_ref, row, col); }
    size_t rows() const
    { return _ref.rows(); }
    size_t cols() const
    { return _ref.cols(); }
private:
    const expr_t &_ref;

    reshaped(const reshaped&);
    //Inaccessible assignment operator.
    reshaped &operator=(const reshaped&);
};

template<typename value_t, typename expr_t, typename tag_t, typename target_shape>
class reshaped<matrix_expression<value_t, expr_t, tag_t>, target_shape, lower_triangular> :
    public base_matrix_expression<
            typename expr_t::value_type,
            reshaped<matrix_expression<value_t, expr_t, tag_t>, target_shape, lower_triangular>,
            tag::expression<util::tag_complexity<tag_t>::value>,
            target_shape
                                  >
{
public:
    typedef typename expr_t::value_type value_type; 
    //typedef typename expr_t::reference reference; ///<Type of reference to element.
    typedef typename expr_t::const_reference const_reference;   

    reshaped(const util::const_ref<expr_t> &a) : _ref(a.ref)
    {   //assertion is needed if user tries to construct reshaped object directly
        ASSERT(reshape_constraint<target_shape>::test(a.ref));
    }
    const_reference operator()(size_t row, size_t col) const
    { return util::translate_access<target_shape, lower_triangular>::get(_ref, row, col); }
    size_t rows() const
    { return _ref.rows(); }
    size_t cols() const
    { return _ref.cols(); }
private:
    const expr_t &_ref;

    reshaped(const reshaped&);
    //Inaccessible assignment operator.
    reshaped &operator=(const reshaped&);
};

namespace util {
template<typename expr_t, typename target_shape, typename side = upper_triangular> struct reshape_return;
template<typename value_t, typename expr_t, typename tag, typename target_shape, typename side>
struct reshape_return<matrix_expression<value_t, expr_t, tag>, target_shape, side>
{
    typedef reshaped<matrix_expression<value_t, expr_t, tag>, target_shape, side> type; 
};
template<typename value_t, typename expr_t, typename tag, typename side>
struct reshape_return<matrix_expression<value_t, expr_t, tag>, rectangular, side>
{
    typedef const expr_t &type; 
};
}//namespace util




template<typename target_shape, typename value_t, typename expr_t, typename tag> inline
typename util::reshape_return<matrix_expression<value_t, expr_t, tag>, target_shape>::type
reshape(const matrix_expression<value_t, expr_t, tag> &a)
{
    reshape_constraint<target_shape>::test(a);
    return util::make_const_ref(get_matrix(a));
}

//reshape<target_shape, side>() form overload
template<typename target_shape, typename side, typename value_t, typename expr_t, typename tag> inline
typename util::reshape_return<matrix_expression<value_t, expr_t, tag>, target_shape, side>::type
reshape(const matrix_expression<value_t, expr_t, tag> &a)
{
    reshape_constraint<target_shape>::test(a);
    return util::make_const_ref(get_matrix(a));
}

namespace util {
template<typename type1, typename type2> struct same_types;
template<typename t>
struct same_types<t, t>
{
    typedef t type;
};
}//namespace util

namespace dummy {
template<typename value_t> inline
const value_t &conj(const value_t &a) { return a; }
template<typename type>
struct operation
{
    typedef type result_type;
    static const result_type &op(const type &a) { return a; }
};
}

namespace op {
template<typename type>
struct plus
{
    typedef type result_type;
    static result_type op(const type &a) { return +a; }
};
template<typename type>
struct minus
{
    typedef type result_type;
    static result_type op(const type &a) { return -a; }
};
template<typename type>
struct conj
{
    typedef type result_type;
    static result_type op(const type &a) { using dummy::conj; return conj(a); }
};
template<typename type1, typename type2>
struct add
{
    typedef typename util::same_types<type1, type2>::type result_type;
    static result_type op(const type1 &a, const type2 &b) { return a + b; }
};
template<typename type1, typename type2>
struct sub
{
    typedef typename util::same_types<type1, type2>::type result_type;
    static result_type op(const type1 &a, const type2 &b) { return a - b; }
};
template<typename type1, typename type2>
struct mul
{
    typedef typename util::same_types<type1, type2>::type result_type;
    static result_type op(const type1 &a, const type2 &b) { return a*b; }
};
template<typename type1, typename type2>
struct div
{
    typedef typename util::same_types<type1, type2>::type result_type;
    static result_type op(const type1 &a, const type2 &b) { return a/b; }
};
}//namespace op


template<typename expression_type, template<typename> class op>
class unary :
    public base_matrix_expression<
        typename expr_traits<expression_type>::type::value_type,
        unary<expression_type, op>,     //by convention this expression itself has constant complexity
        tag::expression<expr_traits<expression_type>::complexity>,
        typename expr_traits<expression_type>::shape
                                  >
{
    typedef typename expr_traits<expression_type>::type expr_t;
public:
    typedef typename expr_t::value_type value_type; 
    typedef const value_type const_reference;   

    unary(const expr_t &a) : _ref(a)
    {}
    const_reference operator()(size_t row, size_t col) const
    { return op<value_type>::op(_ref(row, col)); }
    size_t rows() const
    { return _ref.rows(); }
    size_t cols() const
    { return _ref.cols(); }

private:
    const expr_t &_ref;

    unary(const unary&);    //do not implement
    unary &operator=(const unary&); //do not implement
};




template<typename value_t, typename expr_t, typename tag> inline
const unary<matrix_expression<value_t, expr_t, tag>, op::plus>
operator+(const matrix_expression<value_t, expr_t, tag> &a)
{ return get_matrix(a); }


template<typename value_t, typename expr_t, typename tag, typename shape> inline
const unary<shaped_expression<value_t, expr_t, tag, shape>, op::plus>
operator+(const shaped_expression<value_t, expr_t, tag, shape> &a)
{ return get_matrix(a); }


template<typename value_t, typename expr_t, typename tag> inline
const unary<matrix_expression<value_t, expr_t, tag>, op::minus>
operator-(const matrix_expression<value_t, expr_t, tag> &a)
{ return get_matrix(a); }


template<typename value_t, typename expr_t, typename tag, typename shape> inline
const unary<shaped_expression<value_t, expr_t, tag, shape>, op::minus>
operator-(const shaped_expression<value_t, expr_t, tag, shape> &a)
{ return get_matrix(a); }


template<typename value_t, typename expr_t, typename tag> inline
const unary<matrix_expression<value_t, expr_t, tag>, op::conj>
conj(const matrix_expression<value_t, expr_t, tag> &a)
{ return get_matrix(a); }


template<typename value_t, typename expr_t, typename tag, typename shape> inline
const unary<shaped_expression<value_t, expr_t, tag, shape>, op::conj>
conj(const shaped_expression<value_t, expr_t, tag, shape> &a)
{ return get_matrix(a); }


template<typename expression_type, template<typename> class op = dummy::operation>
class transposed :
    public base_matrix_expression<
        typename expr_traits<expression_type>::type::value_type,
        transposed<expression_type, op>,    //by convention this expression has constant complexity
        tag::expression<expr_traits<expression_type>::complexity>,
        typename expr_traits<expression_type>::shape
                                  >
{
    typedef typename expr_traits<expression_type>::type expr_t;
public:
    typedef typename expr_t::value_type value_type; 
    typedef typename op<value_type>::result_type const_reference;   

    transposed(const expr_t &a) : _ref(a)
    {}
    const_reference operator()(size_t row, size_t col) const
    { return op<value_type>::op(_ref(col, row)); }
    size_t rows() const
    { return _ref.cols(); }
    size_t cols() const
    { return _ref.rows(); }

private:
    const expr_t &_ref;

    transposed(const transposed&);  //do not implement
    transposed &operator=(const transposed&);   //do not implement
};




template<typename value_t, typename expr_t, typename tag> inline
const transposed<matrix_expression<value_t, expr_t, tag> >
transpose(const matrix_expression<value_t, expr_t, tag> &a)
{ return get_matrix(a); }


template<typename value_t, typename expr_t, typename tag, typename shape> inline
const transposed<shaped_expression<value_t, expr_t, tag, shape> >
transpose(const shaped_expression<value_t, expr_t, tag, shape> &a)
{ return get_matrix(a); }


template<typename value_t, typename expr_t, typename tag> inline
const expr_t &transpose(const shaped_expression<value_t, expr_t, tag, symmetric> &a)
{ return get_matrix(a); }


template<typename value_t, typename expr_t, typename tag> inline
const reshaped<shaped_expression<value_t, expr_t, tag, lower_triangular>, upper_triangular, lower_triangular>
transpose(const shaped_expression<value_t, expr_t, tag, lower_triangular> &a)
{ return get_matrix(a); }


template<typename value_t, typename expr_t, typename tag> inline
const reshaped<shaped_expression<value_t, expr_t, tag, upper_triangular>, lower_triangular, upper_triangular>
transpose(const shaped_expression<value_t, expr_t, tag, upper_triangular> &a)
{ return get_matrix(a); }


template<typename expression_type>
class scaled :
    public base_matrix_expression<
        typename expr_traits<expression_type>::type::value_type,
        scaled<expression_type>,
        tag::expression<util::select<expr_traits<expression_type>::complexity, 2>::max>,
        typename expr_traits<expression_type>::shape
                                  >
{
    typedef typename expr_traits<expression_type>::type expr_t;
public:
    typedef typename expr_t::value_type value_type; 
    typedef const value_type const_reference;   

    scaled(const scaled &a) : _ref(a._ref), _scale(a._scale)
    {}
    scaled(const expr_t &a, const_reference s) : _ref(a), _scale(s)
    {}
    const_reference operator()(size_t row, size_t col) const
    { return _ref(row, col)*_scale; }
    size_t rows() const
    { return _ref.rows(); }
    size_t cols() const
    { return _ref.cols(); }

private:
    const expr_t &_ref;
    const value_type _scale;

    scaled &operator=(const scaled&);   //do not implement
};




template<typename value_t, typename expr_t, typename tag> inline
const scaled<matrix_expression<value_t, expr_t, tag> >
operator*(const matrix_expression<value_t, expr_t, tag> &a, const typename expr_t::value_type &s)
{ return scaled<matrix_expression<value_t, expr_t, tag> >(get_matrix(a), s); }


template<typename value_t, typename expr_t, typename tag, typename shape> inline
const scaled<shaped_expression<value_t, expr_t, tag, shape> >
operator*(const shaped_expression<value_t, expr_t, tag, shape> &a, const typename expr_t::value_type &s)
{ return scaled<shaped_expression<value_t, expr_t, tag, shape> >(get_matrix(a), s); }


template<typename value_t, typename expr_t, typename tag> inline
const scaled<matrix_expression<value_t, expr_t, tag> >
operator*(const typename expr_t::value_type &s, const matrix_expression<value_t, expr_t, tag> &a)
{ return scaled<matrix_expression<value_t, expr_t, tag> >(get_matrix(a), s); }


template<typename value_t, typename expr_t, typename tag, typename shape> inline
const scaled<shaped_expression<value_t, expr_t, tag, shape> >
operator*(const typename expr_t::value_type &s, const shaped_expression<value_t, expr_t, tag, shape> &a)
{ return scaled<shaped_expression<value_t, expr_t, tag, shape> >(get_matrix(a), s); }

namespace util {
template<typename shape1, typename shape2>
struct binary_shapes
{
    typedef rectangular resulting_shape;    
};
template<typename shape>
struct binary_shapes<shape, shape>
{
    typedef shape resulting_shape;  
};
template<typename expr1, typename expr2>
struct pair_const_ref
{
    const expr1 &first; 
    const expr2 &second;    

    pair_const_ref(const expr1 &a, const expr2 &b) : first(a), second(b)
    {}

private:
    pair_const_ref &operator=(const pair_const_ref&);   //do not define
};
template<typename expr1, typename expr2>
pair_const_ref<expr1, expr2> make_pair_const_ref(const expr1 &a, const expr2 &b)
{ return pair_const_ref<expr1, expr2>(a, b); }
}//namespace util


template<typename expression1, typename expression2, template<typename, typename> class op>
class binary :
    public base_matrix_expression<
        typename util::same_types<typename expr_traits<expression1>::value_type,
                                  typename expr_traits<expression2>::value_type>::type,
        binary<expression1, expression2, op>,
        tag::expression<util::select<expr_traits<expression1>::complexity,
                                     expr_traits<expression2>::complexity,
                                     2>::max>,
        typename util::binary_shapes<typename expr_traits<expression1>::shape,
                                     typename expr_traits<expression2>::shape>::resulting_shape
                                     >
{
    typedef typename expr_traits<expression1>::type expr1;
    typedef typename expr_traits<expression2>::type expr2;
public:
    typedef typename op<typename expr1::value_type, typename expr2::value_type>::result_type value_type;    
    typedef const value_type const_reference;   

    binary(const util::pair_const_ref<expr1, expr2> &p) :
    _left(p.first), _right(p.second)
    {
        if (_left.rows()!=_right.rows() || _left.cols()!=_right.cols())
            throw size_error("sizes of matrices must be equal");
    }
    const_reference operator()(size_t row, size_t col) const
    { return op<typename expr1::value_type, typename expr2::value_type>::op(_left(row, col), _right(row, col)); }
    size_t rows() const
    {
        ASSERT(_left.rows()==_right.rows());
        return _left.rows();
    }
    size_t cols() const
    {
        ASSERT(_left.cols()==_right.cols());
        return _left.cols();
    }

private:
    const expr1 &_left;
    const expr2 &_right;

    binary(const binary&);  //do not implement
    binary &operator=(const binary&);   //do not implement
};




template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
const binary<matrix_expression<value_t, expr1, tag1>, matrix_expression<value_t, expr2, tag2>, op::add>
operator+(const matrix_expression<value_t, expr1, tag1> &left,
          const matrix_expression<value_t, expr2, tag2> &right
          )
{ return util::make_pair_const_ref(get_matrix(left), get_matrix(right)); }


template<typename value_t, typename expr1, typename tag1, typename shape1, typename expr2, typename tag2, typename shape2> inline
const binary<shaped_expression<value_t, expr1, tag1, shape1>, shaped_expression<value_t, expr2, tag2, shape1>, op::add>
operator+(const shaped_expression<value_t, expr1, tag1, shape1> &left,
          const shaped_expression<value_t, expr2, tag2, shape2> &right
          )
{ return util::make_pair_const_ref(get_matrix(left), get_matrix(right)); }


template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
const binary<matrix_expression<value_t, expr1, tag1>, matrix_expression<value_t, expr2, tag2>, op::sub>
operator-(const matrix_expression<value_t, expr1, tag1> &left,
          const matrix_expression<value_t, expr2, tag2> &right
          )
{ return util::make_pair_const_ref(get_matrix(left), get_matrix(right)); }


template<typename value_t, typename expr1, typename tag1, typename shape1, typename expr2, typename tag2, typename shape2> inline
const binary<shaped_expression<value_t, expr1, tag1, shape1>, shaped_expression<value_t, expr2, tag2, shape2>, op::sub>
operator-(const shaped_expression<value_t, expr1, tag1, shape1> &left,
          const shaped_expression<value_t, expr2, tag2, shape2> &right
          )
{ return util::make_pair_const_ref(get_matrix(left), get_matrix(right)); }


template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
const binary<matrix_expression<value_t, expr1, tag1>, matrix_expression<value_t, expr2, tag2>, op::mul>
mul(const matrix_expression<value_t, expr1, tag1> &left,
    const matrix_expression<value_t, expr2, tag2> &right
    )
{ return util::make_pair_const_ref(get_matrix(left), get_matrix(right)); }


template<typename value_t, typename expr1, typename tag1, typename shape1, typename expr2, typename tag2, typename shape2> inline
const binary<shaped_expression<value_t, expr1, tag1, shape1>, shaped_expression<value_t, expr2, tag2, shape2>, op::mul>
mul(const shaped_expression<value_t, expr1, tag1, shape1> &left,
    const shaped_expression<value_t, expr2, tag2, shape2> &right
    )
{ return util::make_pair_const_ref(get_matrix(left), get_matrix(right)); }


template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
const binary<matrix_expression<value_t, expr1, tag1>, matrix_expression<value_t, expr2, tag2>, op::div>
div(const matrix_expression<value_t, expr1, tag1> &left,
    const matrix_expression<value_t, expr2, tag2> &right
    )
{ return util::make_pair_const_ref(get_matrix(left), get_matrix(right)); }


template<typename value_t, typename expr1, typename tag1, typename shape1, typename expr2, typename tag2, typename shape2> inline
const binary<shaped_expression<value_t, expr1, tag1, shape1>, shaped_expression<value_t, expr2, tag2, shape2>, op::div>
div(const shaped_expression<value_t, expr1, tag1, shape1> &left,
    const shaped_expression<value_t, expr2, tag2, shape2> &right
    )
{ return util::make_pair_const_ref(get_matrix(left), get_matrix(right)); }


template<typename expression_type>
class vector_ref<expression_type, tag::reference, orientation::column> :
    public vector_expression<typename expr_traits<expression_type>::value_type,
                             vector_ref<expression_type, tag::reference, orientation::column>,
                             tag::reference,
                             orientation::column
                             >
{
    typedef typename expr_traits<expression_type>::type expr_t;
public:
    typedef typename expr_t::value_type value_type; 
    typedef typename expr_t::reference reference;   
    typedef typename expr_t::const_reference const_reference;   

    vector_ref(vector_ref &a) : _ref(a._ref), _col(a._col)
    {}
    vector_ref(rvalue<vector_ref> &a) : _ref(a._ref), _col(a._col)
    {}
    vector_ref(expr_t &a, size_t column) : _ref(a), _col(column)
    {}
    reference operator()(size_t i)
    { return _ref(i, _col); }
    const_reference operator()(size_t i) const
    { return _ref(i, _col); }
    reference operator[](size_t i)
    { return _ref(i, _col); }
    const_reference operator[](size_t i) const
    { return _ref(i, _col); }
    size_t size() const
    { return _ref.rows(); }
    //Conversion to "reference to rvalue".
    operator rvalue<vector_ref>&()
    { return to_rvalue_ref(*this); }

private:
    expr_t &_ref;
    const size_t _col;

    vector_ref &operator=(const vector_ref&);   //do not implement

    friend expr_t &get_vector_ref(vector_ref &a) { return a._ref; }
    friend const expr_t &get_vector_ref(const vector_ref &a) { return a._ref; }
    friend size_t get_vector_column(const vector_ref &a) { return a._col; }
};

template<typename expression_type, size_t complexity>
class vector_ref<expression_type, tag::expression<complexity>, orientation::column> :
    public vector_expression<typename expr_traits<expression_type>::value_type,
                             vector_ref<expression_type, tag::expression<complexity>, orientation::column>,
                             tag::expression<util::select<complexity, 1>::max - 1>,
                             orientation::column
                             >
{
    typedef typename expr_traits<expression_type>::type expr_t;
public:
    typedef typename expr_t::value_type value_type; 
    typedef typename expr_t::const_reference const_reference;   

    vector_ref(const vector_ref &a) : _ref(a._ref), _col(a._col)
    {}
    vector_ref(const expr_t &a, size_t col) : _ref(a), _col(col)
    {}
    const_reference operator()(size_t i) const
    { return _ref(i, _col); }
    const_reference operator[](size_t i) const
    { return _ref(i, _col); }
    size_t size() const
    { return _ref.rows(); }

private:
    const expr_t &_ref;
    const size_t _col;

    vector_ref &operator=(const vector_ref&);   //do not implement

    friend const expr_t &get_vector_ref(vector_ref &a) { return a._ref; }
    friend const expr_t &get_vector_ref(const vector_ref &a) { return a._ref; }
    friend size_t get_vector_column(const vector_ref &a) { return a._col; }
};



template<typename expression_type, typename tag> inline
typename util::non_zero_in_column<typename expr_traits<expression_type>::shape>::first_type
first_non_zero(const vector_ref<expression_type, tag, orientation::column> &v)
{ return first_non_zero_in_column(get_vector_ref(v), get_vector_column(v)); }
template<typename expression_type, typename tag> inline
typename util::non_zero_in_column<typename expr_traits<expression_type>::shape>::last_type
last_non_zero(const vector_ref<expression_type, tag, orientation::column> &v)
{ return last_non_zero_in_column(get_vector_ref(v), get_vector_column(v)); }


template<typename expression_type, typename tag>
struct iterator<vector_ref<expression_type, tag, orientation::column> >
{
    typedef typename column_iterator<typename expr_traits<expression_type>::type>::type type;
    static type begin(vector_ref<expression_type, tag, orientation::column> &v, size_t i = 0)
    { return column_iterator<typename expr_traits<expression_type>::type>::begin(get_vector_ref(v), get_vector_column(v), i); }
};

template<typename expression_type, typename tag>
struct const_iterator<vector_ref<expression_type, tag, orientation::column> >
{
    typedef typename column_const_iterator<typename expr_traits<expression_type>::type>::type type;
    static type begin(const vector_ref<expression_type, tag, orientation::column> &v, size_t i = 0)
    { return column_const_iterator<typename expr_traits<expression_type>::type>::begin(get_vector_ref(v), get_vector_column(v), i); }
};

namespace util {
template<typename expression_type>
struct row_ref { typedef vector_ref<expression_type, tag::reference, orientation::row> type; };
template<typename expression_type>
struct const_row_ref { typedef const vector_ref<expression_type, tag::expression<expr_traits<expression_type>::complexity>, orientation::row> type; };
template<typename expression_type>
struct column_ref { typedef vector_ref<expression_type, tag::reference, orientation::column> type; };
template<typename expression_type>
struct const_column_ref { typedef const vector_ref<expression_type, tag::expression<expr_traits<expression_type>::complexity>, orientation::column> type; };
}//namespace util



template<typename value_t, typename expr_t, typename tag> inline
typename util::row_ref<matrix_expression<value_t, expr_t, tag> >::type
row(matrix_expression<value_t, expr_t, tag> &a, size_t r)
{
    ASSERT(r<get_matrix(a).rows());
    return util::row_ref<matrix_expression<value_t, expr_t, tag> >::type(get_matrix(a), r);
}

template<typename value_t, typename expr_t, typename tag_t> inline
typename util::const_row_ref<matrix_expression<value_t, expr_t, tag_t> >::type
row(const matrix_expression<value_t, expr_t, tag_t> &a, size_t r)
{
    ASSERT(r<get_matrix(a).rows());
    return util::const_row_ref<matrix_expression<value_t, expr_t, tag_t> >::type(get_matrix(a), r);
}

template<typename value_t, typename expr_t, typename tag, typename shape> inline
typename util::row_ref<shaped_expression<value_t, expr_t, tag, shape> >::type
row(shaped_expression<value_t, expr_t, tag, shape> &a, size_t r)
{
    ASSERT(r<get_matrix(a).rows());
    return util::row_ref<shaped_expression<value_t, expr_t, tag, shape> >::type(get_matrix(a), r);
}

template<typename value_t, typename expr_t, typename tag_t, typename shape> inline
typename util::const_row_ref<shaped_expression<value_t, expr_t, tag_t, shape> >::type
row(const shaped_expression<value_t, expr_t, tag_t, shape> &a, size_t r)
{
    ASSERT(r<get_matrix(a).rows());
    return util::const_row_ref<shaped_expression<value_t, expr_t, tag_t, shape> >::type(get_matrix(a), r);
}

template<typename value_t, typename expr_t, typename tag> inline
typename util::column_ref<matrix_expression<value_t, expr_t, tag> >::type
column(matrix_expression<value_t, expr_t, tag> &a, size_t c)
{
    ASSERT(c<get_matrix(a).cols());
    return util::column_ref<matrix_expression<value_t, expr_t, tag> >::type(get_matrix(a), c);
}

template<typename value_t, typename expr_t, typename tag_t> inline
typename util::const_column_ref<matrix_expression<value_t, expr_t, tag_t> >::type
column(const matrix_expression<value_t, expr_t, tag_t> &a, size_t c)
{
    ASSERT(c<get_matrix(a).cols());
    return util::const_column_ref<matrix_expression<value_t, expr_t, tag_t> >::type(get_matrix(a), c);
}

template<typename value_t, typename expr_t, typename tag, typename shape> inline
typename util::column_ref<shaped_expression<value_t, expr_t, tag, shape> >::type
column(shaped_expression<value_t, expr_t, tag, shape> &a, size_t c)
{
    ASSERT(c<get_matrix(a).cols());
    return util::column_ref<shaped_expression<value_t, expr_t, tag, shape> >::type(get_matrix(a), c);
}

template<typename value_t, typename expr_t, typename tag_t, typename shape> inline
typename util::const_column_ref<shaped_expression<value_t, expr_t, tag_t, shape> >::type
column(const shaped_expression<value_t, expr_t, tag_t, shape> &a, size_t c)
{
    ASSERT(c<get_matrix(a).cols());
    return util::const_column_ref<shaped_expression<value_t, expr_t, tag_t, shape> >::type(get_matrix(a), c);
}



template<typename expr_t> inline
size_t first_non_zero(const expr_t&)
{ return 0; }
template<typename expr_t> inline
size_t last_non_zero(const expr_t &v)
{ return v.size() - 1; }

namespace util {
struct min { static size_t op(size_t a, size_t b) { return std::min(a, b); } };
struct max { static size_t op(size_t a, size_t b) { return std::max(a, b); } };
template<typename op> inline
size_t combine_indices(size_t a, size_t b)
{ ASSERT(a==b); return a; }
template<typename op> inline
size_t combine_indices(constrained_index a, size_t b)
{ ASSERT(op::op(a, b)==a); return a; }
template<typename op> inline
size_t combine_indices(size_t a, constrained_index b)
{ ASSERT(op::op(a, b)==b); return b; }
template<typename op> inline
size_t combine_indices(constrained_index a, constrained_index b)
{ return op::op(a, b); }
}//namespace util



template<typename expr1, typename expr2>
size_t first_non_zero(const expr1 &a, const expr2 &b)
{ return util::combine_indices<util::max>(first_non_zero(a), first_non_zero(b)); }
template<typename expr1, typename expr2>
size_t last_non_zero(const expr1 &a, const expr2 &b)
{ return util::combine_indices<util::min>(last_non_zero(a), last_non_zero(b)); }





template<template<typename, typename> class op, typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
typename op<value_t, value_t>::result_type
accumulate(const vector_expression<value_t, expr1, tag1> &a, const vector_expression<value_t, expr2, tag2> &b)
{
    typedef typename op<value_t, value_t>::result_type result_type;
    const expr1 &v1 = get_vector(a);
    const expr2 &v2 = get_vector(b);
    ASSERT(v1.size()==v2.size());
    result_type sum = traits<result_type>::zero;
    const size_t size = v1.size();
    /*for (; i<=last; ++i)
        sum = sum + op<value_t, value_t>::op(v1(i), v2(i));*/

    typename const_iterator<expr1>::type i1 = begin(v1, 0);
    typename const_iterator<expr2>::type i2 = begin(v2, 0);
    for (size_t i = 0; i<size; ++i, ++i1, ++i2)
        sum = sum + op<value_t, value_t>::op(*i1, *i2);

    return sum;
}
template<template<typename, typename> class op_t, typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
typename op_t<value_t, value_t>::result_type
accumulate_skip_zeros(const vector_expression<value_t, expr1, tag1> &a, const vector_expression<value_t, expr2, tag2> &b)
{
    typedef op_t<value_t, value_t> op;
    typedef typename op::result_type result_type;
    const expr1 &v1 = get_vector(a);
    const expr2 &v2 = get_vector(b);
    ASSERT(v1.size()==v2.size());
    result_type sum = traits<result_type>::zero;
    size_t i = first_non_zero(v1, v2);
    const size_t last = last_non_zero(v1, v2);
    ASSERT(last<v1.size());
    /*for (; i<=last; ++i)
        sum = sum + op<value_t, value_t>::op(v1(i), v2(i));*/

    typename const_iterator<expr1>::type i1 = begin(v1, i);
    typename const_iterator<expr2>::type i2 = begin(v2, i);
    for (; i<=last; ++i, ++i1, ++i2)
        sum = sum + op::op(*i1, *i2);

    return sum;
}


template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
typename op::mul<value_t, value_t>::result_type
sum_of_products(const vector_expression<value_t, expr1, tag1> &a,
                const vector_expression<value_t, expr2, tag2> &b)
{ return accumulate_skip_zeros<op::mul>(a, b); }

namespace util {
template<typename shape1, typename shape2>
struct product_shapes
{
    typedef rectangular resulting_shape;
};
template<>
struct product_shapes<lower_triangular, lower_triangular>
{
    typedef lower_triangular resulting_shape;
};
template<>
struct product_shapes<upper_triangular, upper_triangular>
{
    typedef upper_triangular resulting_shape;
};
template<typename expression_type, size_t complexity = expr_traits<expression_type>::complexity>
struct const_complexity
{
    typedef const matrix<typename expr_traits<expression_type>::value_type,
                         typename expr_traits<expression_type>::shape> type;
    typedef type &const_reference;
};
template<typename expression_type>
struct const_complexity<expression_type, 0>
{
    typedef const typename expr_traits<expression_type>::type &type;
    typedef type const_reference;
};
}//namespace util


template<typename expression1, typename expression2>
class product :
    public base_matrix_expression<
        typename util::same_types<typename expr_traits<expression1>::value_type,
                                  typename expr_traits<expression2>::value_type>::type,
        product<expression1, expression2>,
        tag::expression<3>,
        typename util::product_shapes<typename expr_traits<expression1>::shape,
                                      typename expr_traits<expression2>::shape>::resulting_shape
                                      >
{
    typedef typename expr_traits<expression1>::type expr1;
    typedef typename expr_traits<expression2>::type expr2;
    typedef typename util::same_types<typename expr1::value_type, typename expr2::value_type>::type deduced_value_type;
public:
    typedef deduced_value_type value_type;  
    typedef const value_type const_reference;   

    product(const util::pair_const_ref<expr1, expr2> &p) :
    _left(p.first), _right(p.second)
    {
        if (_left.cols()!=_right.rows())
            throw size_error("number of columns of the first matrix and number of rows of the second matrix must be equal");
    }

    const_reference operator()(size_t r, size_t c) const
    { return sum_of_products(row(_left, r), column(_right, c)); }
    size_t rows() const
    { return _left.rows(); }
    size_t cols() const
    { return _right.cols(); }
private:
    typename util::const_complexity<expression1>::type _left;
    typename util::const_complexity<expression2>::type _right;

    product(const product&);    //do not implement
    product &operator=(const product&); //do not implement

    friend typename util::const_complexity<expression1>::const_reference left(const product &p) { return p._left; }
    friend typename util::const_complexity<expression2>::const_reference right(const product &p) { return p._right; }
};

template<typename value_t>
template<typename expression1, typename expression2, typename tag>
void matrix<value_t>::assign(const matrix_expression<value_type, product<expression1, expression2>, tag> &rhs)
{
    const product<expression1, expression2> &p = get_matrix(rhs);
    const size_t
        rows = p.rows(),
        cols = p.cols();
    if (rows*cols<64*64 && false)
    {
        this->resize(rows, cols);
        for (size_t i = 0, k = 0; i<_rows; ++i)
            for (size_t j = 0; j<_cols; ++j, ++k)
                _p[k] = p(i, j);
        return;
    }
    typename util::const_complexity<expression1>::const_reference a = left(p);
    typename util::const_complexity<expression2>::const_reference b = right(p);
    typedef typename expr_traits<expression1>::type left_t;
    typedef typename expr_traits<expression2>::type right_t;
    this->resize(a.rows(), b.cols(), traits<value_type>::zero);
    const size_t n = a.cols();
    ASSERT(n==b.rows());
    for (size_t i = 0; i<_rows; ++i)
    {
        const size_t
            k0 = first_non_zero(row(a, i)),
            k1 = last_non_zero(row(a, i));
        ASSERT(k1<n);
        typename row_const_iterator<left_t>::type i1(a, i, k0);
        value_type * const p_i = &_p[i*_cols];
        for (size_t k = k0; k<=k1; ++k, ++i1)
        {
            value_type *p = p_i;
            //const value_type a1 = a(i, k);
            const value_type a1 = *i1;
            const size_t
                first = first_non_zero(row(b, k)),
                last = last_non_zero(row(b, k));
            ASSERT(last<_cols);
            typename row_const_iterator<right_t>::type i2(b, k, first);
            p += first;
            for (size_t j = first; j<=last; ++j, ++p, ++i2)
                *p = *p + a1*(*i2);
        }
    }
    return;
}





template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
product<matrix_expression<value_t, expr1, tag1>, matrix_expression<value_t, expr2, tag2> >
operator*(const matrix_expression<value_t, expr1, tag1> &left,
          const matrix_expression<value_t, expr2, tag2> &right
          )
{ return util::pair_const_ref<expr1, expr2>(get_matrix(left), get_matrix(right)); }


template<typename value_t, typename expr1, typename tag1, typename shape1, typename expr2, typename tag2, typename shape2> inline
product<shaped_expression<value_t, expr1, tag1, shape1>, shaped_expression<value_t, expr2, tag2, shape2> >
operator*(const shaped_expression<value_t, expr1, tag1, shape1> &left,
          const shaped_expression<value_t, expr2, tag2, shape2> &right
          )
{ return util::pair_const_ref<expr1, expr2>(get_matrix(left), get_matrix(right)); }

namespace util {
template<typename value_t, typename expr_t, typename tag_t, size_t complexity>
struct const_complexity<vector_expression<value_t, expr_t, tag_t>, complexity>
{
    typedef const vector<value_t> type;
};
template<typename value_t, typename expr_t, typename tag_t>
struct const_complexity<vector_expression<value_t, expr_t, tag_t>, 0>
{
    typedef const expr_t &type;
};
}//namespace util


template<typename expression1, typename value_t, typename expr2, typename tag2>
class product<expression1, vector_expression<value_t, expr2, tag2> > :
    public vector_expression<
        typename util::same_types<typename expr_traits<expression1>::value_type, value_t>::type,
        product<expression1, vector_expression<value_t, expr2, tag2> >,
        tag::expression<2>,
        typename util::vector_orientation<expr2>::type
                             >
{
    typedef typename expr_traits<expression1>::type expr1;
    typedef typename util::same_types<typename expr1::value_type, typename expr2::value_type>::type deduced_value_type;
public:
    typedef deduced_value_type value_type;  
    typedef const value_type const_reference;   

    product(const util::pair_const_ref<expr1, expr2> &p) :
    _left(p.first), _right(p.second)
    {
        if (_left.cols()!=_right.size())
            throw size_error("number of columns of matrix and size of vector must be equal");
    }

    const_reference operator()(size_t i) const
    { return sum_of_products(row(_left, i), _right); }
    size_t size() const
    { return _left.rows(); }
private:
    typename util::const_complexity<expression1>::type _left;
    typename util::const_complexity<vector_expression<value_t, expr2, tag2> >::type _right;

    product(const product&);    //do not implement
    product &operator=(const product&); //do not implement
};




template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
product<matrix_expression<value_t, expr1, tag1>, vector_expression<value_t, expr2, tag2> >
operator*(const matrix_expression<value_t, expr1, tag1> &left,
          const vector_expression<value_t, expr2, tag2, orientation::column> &right
          )
{ return util::pair_const_ref<expr1, expr2>(get_matrix(left), get_vector(right)); }


template<typename value_t, typename expr1, typename tag1, typename shape, typename expr2, typename tag2> inline
product<shaped_expression<value_t, expr1, tag1, shape>, vector_expression<value_t, expr2, tag2> >
operator*(const shaped_expression<value_t, expr1, tag1, shape> &left,
          const vector_expression<value_t, expr2, tag2, orientation::column> &right
          )
{ return util::pair_const_ref<expr1, expr2>(get_matrix(left), get_vector(right)); }

template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
product<matrix_expression<value_t, expr1, tag1>, vector_expression<value_t, expr2, tag2> >
map(const matrix_expression<value_t, expr1, tag1> &linear_map,
    const vector_expression<value_t, expr2, tag2> &v
    )
{ return util::pair_const_ref<expr1, expr2>(get_matrix(linear_map), get_vector(v)); }


template<typename value_t, typename expr1, typename tag1, typename shape, typename expr2, typename tag2> inline
product<shaped_expression<value_t, expr1, tag1, shape>, vector_expression<value_t, expr2, tag2> >
map(const shaped_expression<value_t, expr1, tag1, shape> &linear_map,
    const vector_expression<value_t, expr2, tag2> &v
    )
{ return util::pair_const_ref<expr1, expr2>(get_matrix(linear_map), get_vector(v)); }


template<typename value_t, typename expr1, typename expr2, typename tag2> inline
expr1 &operator+=(matrix_expression<value_t, expr1, tag::storage> &left,
                  const matrix_expression<value_t, expr2, tag2> &right
                  )
{
    expr1 &a = get_matrix(left);
    expr1(a + right).swap(a);
    return a;
}


template<typename value_t, typename expr1, typename shape, typename expr2, typename tag2> inline
expr1 &operator+=(shaped_expression<value_t, expr1, tag::storage, shape> &left,
                  const shaped_expression<value_t, expr2, tag2, shape> &right
                  )
{
    expr1 &a = get_matrix(left);
    expr1(a + right).swap(a);
    return a;
}


template<typename value_t, typename expr1, typename expr2, typename tag2> inline
expr1 &operator-=(matrix_expression<value_t, expr1, tag::storage> &left,
                  const matrix_expression<value_t, expr2, tag2> &right
                  )
{
    expr1 &a = get_matrix(left);
    expr1(a - right).swap(a);
    return a;
}


template<typename value_t, typename expr1, typename shape, typename expr2, typename tag2> inline
expr1 &operator-=(shaped_expression<value_t, expr1, tag::storage, shape> &left,
                  const shaped_expression<value_t, expr2, tag2, shape> &right
                  )
{
    expr1 &a = get_matrix(left);
    expr1(a - right).swap(a);
    return a;
}


template<typename value_t, typename expr_t> inline
expr_t &operator*=(matrix_expression<value_t, expr_t, tag::storage> &left,
                  const value_t &right
                  )
{
    expr_t &a = get_matrix(left);
    const size_t
        rows = a.rows(),
        cols = a.cols();
    for (size_t i = 0; i<rows; ++i)
        for (size_t j = 0; j<cols; ++j)
            a(i, j) *= right;
    return a;
}


template<typename value_t, typename expr_t, typename shape> inline
expr_t &operator*=(shaped_expression<value_t, expr_t, tag::storage, shape> &left,
                  const value_t &right
                  )
{
    expr_t &a = get_matrix(left);
    expr_t(a*right).swap(a);
    return a;
}


template<typename value_t, typename expr1, typename expr2, typename tag2> inline
expr1 &operator*=(matrix_expression<value_t, expr1, tag::storage> &left,
                  const matrix_expression<value_t, expr2, tag2> &right
                  )
{
    expr1 &a = get_matrix(left);
    expr1(a*right).swap(a);
    return a;
}


template<typename value_t, typename expr1, typename expr2, typename tag2> inline
expr1 &operator*=(shaped_expression<value_t, expr1, tag::storage, lower_triangular> &left,
                  const shaped_expression<value_t, expr2, tag2, lower_triangular> &right
                  )
{
    expr1 &a = get_matrix(left);
    expr1(a*right).swap(a);
    return a;
}


template<typename value_t, typename expr1, typename expr2, typename tag2> inline
expr1 &operator*=(shaped_expression<value_t, expr1, tag::storage, upper_triangular> &left,
                  const shaped_expression<value_t, expr2, tag2, upper_triangular> &right
                  )
{
    expr1 &a = get_matrix(left);
    expr1(a*right).swap(a);
    return a;
}




template<typename value_t, typename expr_t> inline
void swap(vector_expression<value_t, expr_t, tag::storage> &a,
          vector_expression<value_t, expr_t, tag::storage> &b)
{ get_vector(a).swap(get_vector(b)); return; }


template<typename value_t, typename expr1, typename expr2> inline
void swap(vector_expression<value_t, expr1, tag::storage> &a,
          vector_expression<value_t, expr2, tag::storage> &b)
{
    expr1 &v1 = get_vector(a);
    expr2 &v2 = get_vector(b);
    expr1 temp = v2;
    v2 = v1;
    swap(v1, temp);
    return;
}

namespace util {
template<typename expr1, typename expr2> inline
void swap_vectors(expr1& a, expr2 &b)
{
    using std::swap;
    const size_t count = std::min(a.size(), b.size());
    for (size_t i = 0; i<count; ++i)
        swap(a(i), b(i));
    return;
}
}//namespace util




template<typename value_t, typename expr1, typename expr2, typename orientation> inline
void swap(vector_expression<value_t, expr1, tag::storage> &a,
          vector_ref<expr2, tag::reference, orientation> b)
{   //ensuring the value type is the same
    typedef typename util::same_types<value_t, expr_traits<expr2>::value_type>::type value_type;
    util::swap_vectors(get_vector(a), b);
    return;
}

template<typename expr1, typename orientation, typename value_t, typename expr2> inline
void swap(vector_ref<expr1, tag::reference, orientation> a,
          vector_expression<value_t, expr2, tag::storage> &b)
{   //ensuring the value type is the same
    typedef typename util::same_types<value_t, expr_traits<expr1>::value_type>::type value_type;
    util::swap_vectors(a, get_vector(b));
    return;
}

template<typename expr1, typename orientation1, typename expr2, typename orientation2> inline
void swap(vector_ref<expr1, tag::reference, orientation1> a,
          vector_ref<expr2, tag::reference, orientation2> b)
{   //ensuring the value type is the same
    typedef typename util::same_types<expr_traits<expr1>::value_type, expr_traits<expr2>::value_type>::type value_type;
    util::swap_vectors(a, b);
    return;
}

namespace util {
template<typename type1, typename type2>
struct mul_conj
{
    typedef typename util::same_types<type1, type2>::type result_type;
    static result_type op(const type1 &a, const type2 &b) { using dummy::conj; return a*conj(b); }
};
}//namespace util


template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
typename util::mul_conj<value_t, value_t>::result_type
dot_product(const vector_expression<value_t, expr1, tag1> &v1,
            const vector_expression<value_t, expr2, tag2> &v2)
{ return accumulate_skip_zeros<util::mul_conj>(v1, v2); }



template<template<typename> class op, typename value_t, typename expr_t, typename tag> inline
typename op<value_t>::result_type
accumulate(const vector_expression<value_t, expr_t, tag> &a)
{
    typedef typename op<value_t>::result_type result_type;
    const expr_t v = get_vector(a);
    result_type sum = traits<result_type>::zero;
    const size_t size = v.size();
    for (size_t i = 0; i<size; ++i)
        sum = sum + op<value_t>::op(v(i));
    return sum;
}
template<template<typename> class op_t, typename value_t, typename expr_t, typename tag> inline
typename op_t<value_t>::result_type
accumulate_skip_zeros(const vector_expression<value_t, expr_t, tag> &a)
{
    //typedef typename op<value_t>::result_type result_type;
    //const expr_t v = get_vector(a);
    //result_type sum = traits<result_type>::zero;
    //const size_t size = v.size();
    //for (size_t i = 0; i<size; ++i)
    //  sum = sum + op<value_t>::op(v(i));
    //return sum;

    typedef op_t<value_t> op;
    typedef typename op::result_type result_type;
    const expr_t &v = get_vector(a);
    result_type sum = traits<result_type>::zero;
    size_t i = first_non_zero(v);
    const size_t last = last_non_zero(v);
    ASSERT(last<v.size());

    typename const_iterator<expr_t>::type it = begin(v, i);
    for (; i<=last; ++i, ++it)
        sum = sum + op::op(*it);

    return sum;
}


namespace util {
template<typename type>
struct sqr
{
    typedef type result_type;
    static result_type op(const type &a) { using dummy::conj; return a*conj(a); }
};
}//namespace util




template<typename value_t, typename expr_t, typename tag>
inline value_t sqrnorm(const vector_expression<value_t, expr_t, tag> &v)
{ return accumulate_skip_zeros<util::sqr>(v); }


template<typename value_t, typename expr_t, typename tag>
inline value_t norm(const vector_expression<value_t, expr_t, tag> &v)
{ return sqrt(sqrnorm(v)); }


template<typename value_t, typename expr_t, typename tag_t, template<typename> class op>
class unary<vector_expression<value_t, expr_t, tag_t>, op> :
    public vector_expression<value_t,
                             unary<vector_expression<value_t, expr_t, tag_t>, op>,
                             tag::expression< util::select<util::tag_complexity<tag_t>::value, 0>::max>,
                             typename util::vector_orientation<expr_t>::type
                             >
{
public:
    typedef value_t value_type; 
    typedef const value_type const_reference;   

    unary(const expr_t &a) : _ref(a)
    {}
    const_reference operator()(size_t i) const
    { return op<value_type>::op(_ref(i)); }
    size_t size() const
    { return _ref.size(); }

private:
    const expr_t &_ref;

    unary(const unary&);    //do not implement
    unary &operator=(const unary&); //do not implement
};




template<typename value_t, typename expr_t, typename tag> inline
const unary<vector_expression<value_t, expr_t, tag>, op::plus>
operator+(const vector_expression<value_t, expr_t, tag> &a)
{ return get_vector(a); }


template<typename value_t, typename expr_t, typename tag> inline
const unary<vector_expression<value_t, expr_t, tag>, op::minus>
operator-(const vector_expression<value_t, expr_t, tag> &a)
{ return get_vector(a); }


template<typename value_t, typename expr_t, typename tag> inline
const unary<vector_expression<value_t, expr_t, tag>, op::conj>
conj(const vector_expression<value_t, expr_t, tag> &a)
{ return get_vector(a); }


template<typename value_t, typename expr_t, typename tag_t>
class scaled<vector_expression<value_t, expr_t, tag_t> > :
    public vector_expression<value_t,
                             scaled<vector_expression<value_t, expr_t, tag_t> >,
                             tag::expression<util::select<util::tag_complexity<tag_t>::value, 2>::max>,
                             typename util::vector_orientation<expr_t>::type
                             >
{
public:
    typedef value_t value_type; 
    typedef const value_type const_reference;   

    scaled(const scaled &a) : _ref(a._ref), _scale(a._scale)
    {}
    scaled(const expr_t &a, const_reference s) : _ref(a), _scale(s)
    {}
    const_reference operator()(size_t i) const
    { return _ref(i)*_scale; }
    size_t size() const
    { return _ref.size(); }

private:
    const expr_t &_ref;
    const value_type _scale;

    scaled &operator=(const scaled&);   //do not implement
};




template<typename value_t, typename expr_t, typename tag> inline
const scaled<vector_expression<value_t, expr_t, tag> >
operator*(const vector_expression<value_t, expr_t, tag> &a, const typename expr_t::value_type &s)
{ return scaled<vector_expression<value_t, expr_t, tag> >(get_vector(a), s); }


template<typename value_t, typename expr_t, typename tag> inline
const scaled<vector_expression<value_t, expr_t, tag> >
operator*(const typename expr_t::value_type &s, const vector_expression<value_t, expr_t, tag> &a)
{ return scaled<vector_expression<value_t, expr_t, tag> >(get_vector(a), s); }

namespace util {
template<typename orienation1, typename orienation2>
struct orientations
{
    typedef orientation::undefined resulting_orientation;
};
template<typename orienation_t>
struct orientations<orienation_t, orienation_t>
{
    typedef orienation_t resulting_orientation;
};
template<typename orienation_t>
struct orientations<orientation::undefined, orienation_t>
{
    typedef orienation_t resulting_orientation;
};
template<typename orienation_t>
struct orientations<orienation_t, orientation::undefined>
{
    typedef orienation_t resulting_orientation;
};
}//namespace util


template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2, template<typename, typename> class op>
class binary<vector_expression<value_t, expr1, tag1>, vector_expression<value_t, expr2, tag2>, op> :
    public vector_expression<
        value_t,
        binary<vector_expression<value_t, expr1, tag1>, vector_expression<value_t, expr2, tag2>, op>,
        tag::expression<util::select<util::tag_complexity<tag1>::value,
                                     util::tag_complexity<tag2>::value,
                                     1>::max>,
        typename util::orientations<typename util::vector_orientation<expr1>::type,
                                    typename util::vector_orientation<expr2>::type>::resulting_orientation
                                  >
{
public:
    typedef typename op<typename expr1::value_type, typename expr2::value_type>::result_type value_type;    
    typedef const value_type const_reference;   

    binary(const util::pair_const_ref<expr1, expr2> &p) :
    _left(p.first), _right(p.second)
    {
        if (_left.size()!=_right.size())
            throw size_error("sizes of vectors must be equal");
    }
    const_reference operator()(size_t i) const
    { return op<typename expr1::value_type, typename expr2::value_type>::op(_left(i), _right(i)); }
    size_t size() const
    {
        ASSERT(_left.size()==_right.size());
        return _left.size();
    }

private:
    const expr1 &_left;
    const expr2 &_right;

    binary(const binary&);  //do not implement
    binary &operator=(const binary&);   //do not implement
};




template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
const binary<vector_expression<value_t, expr1, tag1>, vector_expression<value_t, expr2, tag2>, op::add>
operator+(const vector_expression<value_t, expr1, tag1> &left,
          const vector_expression<value_t, expr2, tag2> &right
          )
{ return util::pair_const_ref<expr1, expr2>(get_vector(left), get_vector(right)); }


template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
const binary<vector_expression<value_t, expr1, tag1>, vector_expression<value_t, expr2, tag2>, op::sub>
operator-(const vector_expression<value_t, expr1, tag1> &left,
          const vector_expression<value_t, expr2, tag2> &right
          )
{ return util::pair_const_ref<expr1, expr2>(get_vector(left), get_vector(right)); }


template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
const binary<vector_expression<value_t, expr1, tag1>, vector_expression<value_t, expr2, tag2>, op::mul>
mul(const vector_expression<value_t, expr1, tag1> &left,
    const vector_expression<value_t, expr2, tag2> &right
    )
{ return util::pair_const_ref<expr1, expr2>(get_vector(left), get_vector(right)); }


template<typename value_t, typename expr1, typename tag1, typename expr2, typename tag2> inline
const binary<vector_expression<value_t, expr1, tag1>, vector_expression<value_t, expr2, tag2>, op::div>
div(const vector_expression<value_t, expr1, tag1> &left,
    const vector_expression<value_t, expr2, tag2> &right
    )
{ return util::pair_const_ref<expr1, expr2>(get_vector(left), get_vector(right)); }


template<typename value_t, typename expr1, typename expr2, typename tag2> inline
expr1 &operator+=(vector_expression<value_t, expr1, tag::storage> &left,
                  const vector_expression<value_t, expr2, tag2> &right
                  )
{
    expr1 &a = get_vector(left);
    expr1(a + get_vector(right)).swap(a);
    return a;
}


template<typename value_t, typename expr1, typename expr2, typename tag2> inline
expr1 &operator+=(vector_expression<value_t, expr1, tag::reference> &left,
                  const vector_expression<value_t, expr2, tag2> &right
                  )
{
    expr1 &a = get_vector(left);
    const expr2 &b = get_vector(right);
    const size_t size = a.size();
    ASSERT(size==b.size());
    for (size_t i = 0; i<size; ++i)
        a(i) += b(i);
    return a;
}


template<typename value_t, typename expr1, typename expr2, typename tag2> inline
expr1 &operator-=(vector_expression<value_t, expr1, tag::storage> &left,
                  const vector_expression<value_t, expr2, tag2> &right
                  )
{
    expr1 &a = get_vector(left);
    expr1(a - get_vector(right)).swap(a);
    return a;
}


template<typename value_t, typename expr1, typename expr2, typename tag2> inline
expr1 &operator-=(vector_expression<value_t, expr1, tag::reference> &left,
                  const vector_expression<value_t, expr2, tag2> &right
                  )
{
    expr1 &a = get_vector(left);
    const expr2 &b = get_vector(right);
    const size_t size = a.size();
    ASSERT(size==b.size());
    for (size_t i = 0; i<size; ++i)
        a(i) -= b(i);
    return a;
}

namespace util {
template<typename type> inline
type &scale_vector(type &a, typename const type::value_type &b)
{
    const size_t size = a.size();
    for (size_t i = 0; i<size; ++i)
        a(i) = a(i)*b;
    return a;
}
}//namespace util




template<typename value_t, typename expr_t> inline
expr_t &operator*=(vector_expression<value_t, expr_t, tag::storage> &left,
                  const value_t &right
                  )
{ return util::scale_vector(get_vector(left), right); }


template<typename value_t, typename expr_t> inline
expr_t &operator*=(vector_expression<value_t, expr_t, tag::reference> &left,
                  const value_t &right
                  )
{ return util::scale_vector(get_vector(left), right); }

}//namespace nabla
#pragma warning(pop)
#endif//TOOLS_MATRIX_INCLUDED