vector_type solve(const matrix_type& A, const vector_type& y)
    {
        typedef typename matrix_type::size_type size_type;
        typedef typename matrix_type::value_type value_type;

        namespace ublas = boost::numeric::ublas;

        matrix_type Q(A.size1(), A.size2()), R(A.size1(), A.size2());

        qr (A, Q, R);

        vector_type b = prod(trans(Q), y);

        vector_type result;
        if (R.size1() > R.size2())
        {
            size_type min = (R.size1() < R.size2() ? R.size1() : R.size2());

            result = ublas::solve(subrange(R, 0, min, 0, min),
                                  subrange(b, 0, min),
                                  ublas::upper_tag());
        }
        else
        {
            result = ublas::solve(R, b, ublas::upper_tag());
        }
        return result;
    }

  /** 
   * 13. 
   * @brief Copy the contents of a matrix.
   * @param[out] A A is overwritten with B.
   * @param[in] B The matrix to be copied.
   */
  static void copy (matrix_type &A, const matrix_type &B) { 

    /** Result matrices should always be of the right size */
    A.Resize (B.Height(), B.Width());

    /** This one is pretty simple, but the order is different */
    elem::Copy (B, A);
  }

  /** 
   * 15. 
   * @brief Extract the diagonal elements of a matrix.
   * @param[in] A The (square) matrix whose diagonal is to be extracted.
   * @param[out] B A diagonal matrix containing entries from A.
   */
  static void diag(const matrix_type& A,
                   matrix_type& B) {
    /* create a zeros matrix of the right dimension and assign to B */
    const int n = A.Width();
    B = zeros(n, n);

    /* Set the diagonals of B */
    for (int i=0; i<n; ++i) B.Set(i, i, A.Get(i,i));
  }

  /** 
   * 6. 
   * @brief Compute the matrix transpose.
   * @param[in] A The matrix to be transposed.
   * @param[out] B B is overwritten with A^{T}
   */
  static void transpose (const matrix_type& A, 
                         matrix_type& B) {

    /** Result matrices should always have sufficient space */
    B.Resize(A.Width(), A.Height());

    /** Compute transpose */
    elem::Transpose (A, B);
  }

  /** 
   * 5. 
   * @brief Multiply one matrix with another.
   * @param[in] A The first matrix
   * @param[in] B The second matrix
   * @param[out] C C is overwritten with (A*B)
   */
  static void multiply (const matrix_type& A,
                        const matrix_type& B,
                        matrix_type& C) {

    /** Result matrices should always have sufficient space */
    C.Resize(A.Height(), B.Width());

    /** We have to do a Gemm */
    elem::Gemm (elem::NORMAL, elem::NORMAL, 1.0, A, B, 0.0, C);
  }

	/** @brief Constructs a tensor with a matrix
	 *
	 * \note Initially the tensor will be two-dimensional.
	 *
	 *  @param v matrix to be copied.
	 */
	BOOST_UBLAS_INLINE
	tensor (const matrix_type &v)
		: tensor_expression_type<self_type>()
		, extents_ ()
		, strides_ ()
		, data_    (v.data())
	{
		if(!data_.empty()){
			extents_ = extents_type{v.size1(),v.size2()};
			strides_ = strides_type(extents_);
		}
	}

// Vectorizes (stacks columns of) an input matrix.
state_type Vec(const matrix_type &A)
{
    state_type vectorized(A.size1()*A.size2());
    
    for (int n = 0; n < A.size1(); n++) {
        for (int m = 0; m < A.size2(); m++) {
            vectorized(n*A.size1()+m) = A(m,n);
        }
    }
    
    return vectorized;
}

  /** 
   * 7. 
   * @brief Compute element-wise negation of a matrix.
   * @param[in] A The matrix to be negated.
   * @param[out] B B is overwritten with -1.0*A
   */
  static void negation (const matrix_type& A,
                        matrix_type& B) {

    /** Result matrices should always have sufficient space */
    B.Resize(A.Height(), A.Width());

    /** Copy over the matrix */
    elem::Copy (A, B);

    /** Multiply by -1.0 */
    elem::Scal(-1.0, B);
  }

	static void calc(matrix_type& a, vector_type& v, const char& jobz) {
		namespace impl = boost::numeric::bindings::lapack::detail;
		//impl::syev(jobz,a,v);

		char      uplo  = 'L';
		integer_t n     = a.size2();
		integer_t lda   = a.size1();
		integer_t lwork =-1;
		integer_t info  = 0;
		value_type dlwork;
		impl::syev(jobz, uplo, n, mtraits::data(a), lda, vtraits::data(v), &dlwork, lwork, info);
	}

    /**
     * Apply columnwise the sketching transform that is described by the
     * the transform with output sketch_of_A.
     */
    void apply (const matrix_type& A,
                output_matrix_type& sketch_of_A,
                columnwise_tag dimension) const {

        const value_type *a = A.LockedBuffer();
        El::Int lda = A.LDim();
        value_type *sa = sketch_of_A.Buffer();
        El::Int ldsa = sketch_of_A.LDim();

        for (El::Int j = 0; j < A.Width(); j++)
            for (El::Int i = 0; i < data_type::_S; i++)
                sa[j * ldsa + i] = a[j * lda + data_type::_samples[i]];
    }

  /** 
   * 4. 
   * @brief Subtract one matrix from another.
   * @param[in] A The first matrix
   * @param[in] B The second matrix
   * @param[out] C C is overwritten with (A-B)
   */
  static void minus (const matrix_type& A,
                     const matrix_type& B, 
                     matrix_type& C) {

    /** Result matrices should always have sufficient space */
    C.Resize(A.Height(), A.Width());

    /** first copy the matrix over */
    elem::Copy (A, C);

    /** now, subtract the other matrix in */
    elem::Axpy (-1.0, B, C);
  }

    bool is_positive_definite( const matrix<T,N,A>& m )
    {
        typedef matrix<T,N,A> matrix_type;
        typedef typename matrix_type::range_type range_type;

        if ( m.row() != m.col() ) return false;
        
        for ( std::size_t i = 1; i != m.row(); ++i )
        {
            const matrix_type a{ m, range_type{0,i}, range_type{0,i} };
            if ( a.det() <= T(0) ) return false;
        }

       return true; 
    }

        // --------------------------------------------------------------------
        static matrix_type _compute_rooted_fa(const matrix_type & fa)
        {
            const size_t K = fa.get_height();
            const size_t J = fa.get_width();

            assert(K > 1);

            matrix_type rooted_fa (K - 1, J);

            for (size_t k = 0; k + 1 < K; k++)
                for (size_t j = 0; j < J; j++)
                    rooted_fa(k, j) = fa(k + 1, j) - fa(0, j);

            return rooted_fa;
        }

    /**
     * Apply the sketching transform on A and write to sketch_of_A.
     * Implementation for columnwise.
     */
    void apply_impl_vdist(const matrix_type& A,
        output_matrix_type& sketch_of_A,
        skylark::sketch::columnwise_tag tag) const {

        // Just a local operation on the Matrix
        _local.apply(A.LockedMatrix(), sketch_of_A.Matrix(), tag);
    }

// Tensor product. Probably a better way of doing this.
matrix_type Kronecker(const matrix_type &A, const matrix_type &B)
{
    matrix_type C(A.size1()*B.size1(), A.size2()*B.size2());
    
    for (int i=0; i < A.size1(); i++) {
        for (int j=0; j < A.size2(); j++) {
            for (int k=0; k < B.size1(); k++) {
                for (int l=0; l < B.size2(); l++) {
                    C(i*B.size1()+k, j*B.size2()+l) = A(i,j)*B(k,l);
                }
            }
        }
    }
    
    return C;
}

	static scalar_type r(const matrix_type& m)
	{
		BOOST_STATIC_ASSERT(Row >= 0);
		BOOST_STATIC_ASSERT(Row < rows);
		BOOST_STATIC_ASSERT(Col >= 0);
		BOOST_STATIC_ASSERT(Col < cols);
		return m.at(Col, Row);
	}

	static scalar_type ir(int row, int col, const matrix_type& m)
	{
		BOOST_ASSERT(row >= 0);
		BOOST_ASSERT(row < rows);
		BOOST_ASSERT(col >= 0);
		BOOST_ASSERT(col < cols);
		return m.at(col, row);
	}

void fill_from_file(matrix_type& tens, const std::string& str)
{
   std::ifstream ifs(str);
   for(int j = 0 ; j < tens.extent<1>(); ++j)
      for(int i = 0 ; i < tens.extent<0>(); ++i)
   {
      ifs >> tens.at(i,j);
   }
}

            // ----------------------------------------------------------------
            explicit q_data(const matrix_type & d)
                : d_ij ()
                , d_ik ()
                , d_jk ()
                , i    ()
                , j    ()
            {
                static const auto k_0_5 = value_type(0.5);

                const auto n     = d.get_height();
                const auto k_n_2 = value_type(n - 2);

                //
                // Cache the row sums; column sums would work equally well
                // because the matrix is symmetric.
                //
                std::vector<value_type> sigma;
                for (size_t c = 0; c < n; c++)
                    sigma.push_back(d.get_row_sum(c));

                //
                // Compute the values of the Q matrix.
                //
                matrix_type q (n, n);
                for (size_t r = 0; r < n; r++)
                    for (size_t c = 0; c < r; c++)
                        q(r, c) = k_n_2 * d(r, c) - sigma[r] - sigma[c];

                //
                // Find the cell with the minimum value.
                //
                i = 1, j = 0;
                for (size_t r = 2; r < n; r++)
                    for (size_t c = 0; c < r; c++)
                        if (q(r, c) < q(i, j))
                            i = r, j = c;

                //
                // Compute distances between the new nodes.
                //
                d_ij = d(i, j);
                d_ik = k_0_5 * (d_ij + ((sigma[i] - sigma[j]) / k_n_2));
                d_jk = d_ij - d_ik;
            }

box<dimension>::box(matrix_type const& edges)
{
    if (edges.size1() != dimension || edges.size2() != dimension) {
        throw std::invalid_argument("edge vectors have invalid dimensionality");
    }
    edges_ = zero_matrix_type(dimension, dimension);
    for (unsigned int i = 0; i < dimension; ++i) {
        edges_(i, i) = edges(i, i);
    }
    if (norm_inf(edges_ - edges) != 0) {
        throw std::invalid_argument("non-cuboid geomtries are not implemented");
    }
    for (unsigned int i = 0; i < dimension; ++i) {
        length_[i] = edges_(i, i);
    }
    length_half_ = 0.5 * length_;

    LOG("edge lengths of simulation domain: " << length_);
}