matrix_vector_operation.h

Go to the documentation of this file.

#ifndef MATRIX_OPERATION_H
#define MATRIX_OPERATION_H

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>

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

#include <stdexcept>

namespace matt_boost {
    namespace ublas = boost::numeric::ublas;

    template<class V>
        BOOST_UBLAS_INLINE
        V trace(ublas::matrix<V>& e) {
            V tr = 0;
            if (e.size1() == e.size2()) {
                ublas::matrix_vector_range<ublas::matrix<V> > diag(e,ublas::range(0,e.size1()),ublas::range(0,e.size2()));
                tr = sum(diag);
            }
            else
                throw std::length_error("Error in function trace() of matrix_vector_operation.h: Wrong matrix dimensions.");
            
            return tr;
        }

    template<class V>
        BOOST_UBLAS_INLINE
        ublas::vector<V> cross_prod(ublas::vector<V>& v1, ublas::vector<V>& v2) {
            ublas::vector<V> v(v1.size());
            if (v1.size() == v2.size() && v1.size() == 3) {
                v(0) = v1(1) * v2(2) - v1(2) * v2(1);
                v(1) = v1(2) * v2(0) - v1(0) * v2(2);
                v(2) = v1(0) * v2(1) - v1(1) * v2(0);
            }
            else
                throw std::length_error("Error in function cross_prod() of matrix_vector_operation.h: Wrong vector dimensions.");
            
            return v;
        }

    template<class V>
        BOOST_UBLAS_INLINE
        ublas::matrix<V> taylor_exp(const ublas::matrix<V>& F, const V ds, const unsigned int order) {
            double fac = 1.0;
            ublas::matrix<V> Fn = ublas::identity_matrix<V>(6);
            ublas::matrix<V> M = Fn;

            for (unsigned int k = 1; k < order; ++k) {
                fac *= ds / V(k);
                Fn = prod(Fn,F);
                M = M + fac * Fn;
            }
            return M;
        }

    template<class M, class E1, class E2, class E3>
        BOOST_UBLAS_INLINE
        M gemmm(const ublas::matrix_expression<E1>& e1, const ublas::matrix_expression<E2>& e2, const ublas::matrix_expression<E3>& e3) {
            M tmp = prod(e2,e3);
            return prod(e1,tmp);
        }
}

#endif