#include <boost/numeric/ublas/hermitian.hpp>

using namespace boost::numeric::ublas;

template< class ValueT >
class four_indexed_container
{
public:
    typedef ValueT value_type;
    typedef hermitian_matrix<value_type, upper, column_major> matrix_type;
    typedef typename matrix_type::size_type size_type;
    typedef typename matrix_type::reference reference;
    typedef typename matrix_type::const_reference const_reference;

    four_indexed_container (const size_type n)
        : size_(n),
          data1_(matrix_size(size_), matrix_size(size_)),
          data2_(matrix_size(size_), matrix_size(size_))
    {}

    ~four_indexed_container()
    {}

    //--------------------------------------------------------------------------
    //  Index access
    //--------------------------------------------------------------------------
    reference operator() (const size_type i, const size_type j,
                          const size_type k, const size_type l)
    {
        assert (i < size_);
        assert (j < size_);
        assert (k < size_);
        assert (l < size_);
        assert (i != j && k != l);
        // m(i,j|k,l) = m(j,i|l,k)
        return i < j ? get_ref(i, j, k, l) : get_ref(j, i, l, k);
    }

private:
    // Helper functions
    static
    size_type matrix_size (size_type n)
    {
        return n * (n - 1) / 2;
    }

    // reference
    reference get_ref(size_type i, size_type j,
                      size_type k, size_type l)
    {
        assert(i < j && k != l);
        size_type ij = i + j * (j - 1) / 2;
        if (k < l)
        {
            size_type kl = k + l * (l - 1) / 2;
            return data1_(ij, kl);
        }
        else
        {
            size_type kl = l + k * (k - 1) / 2;
            return data2_(ij, kl);
        }
    }

private:
    size_type    size_;
    matrix_type  data1_;        // (i,j|k,l) elements with i < j and k < l.
    matrix_type  data2_;        // (i,j|k,l) elements with i < j and k > l.
};



int main()
{
    four_indexed_container<std::complex<double> > a(5);
    // next line should modify a.data1_(0, 0).
    a(0, 1, 0, 1) = std::complex<double>(0, 0); // <-- assertion fails!!

    return 0;
}