regarding element-by-element operations in boost::ublas

Question

I find that boost::ublas doesn't support element-by-element operations and operations in sequence very well (but the efficiency is pretty high :)) I am trying to

D = A^2 .* B^3 .* C

where A, B, C are all square matrices of the same size, operator ".*" denotes the element-by-element operation and ^ is the power of the matrix. With boost:ublas, I wrote

for (int n=0; n<300; n++)
{
  for (int k=0; k<300; k++)
  {
    D(n, k) = pow(abs(A(n, k)), 2)*pow(abs(B(n, k)), 3)*C(n, k);
  }
}

In my program I have many sequent operations like those shown above, anyway I can get the same result but using one line of code instead of loop?

Also, I observe that it seems no valid to assign a constant to all elements of matrix or vector like

boost::numeric::ublas::vector v(100); v = 0.2;

Instead, I have to use a loop to do the assignment again, any better way to save some code? My algorithm is a really long and there are so many tedious operations like those stated above. I tried another numerical library Armadillo, which provides a good way to simply the operations but it currently doesn't suppose sparse matrices (which will spend about 10 times in running my code).

Answer 1

You can assign constant to a vector or matrix easily enough:

vector<double> v = scalar_vector<double>(100, 0.2);

regular vector (but not c_vector or bounded_vector) even has a constructor:

vector<double> v(100, 0.2);

As for element operations, you can also easily define your own, it' just the boilerplate code, e.g. for absolute power:) :

namespace boost { namespace numeric { namespace ublas {
template<class M, class P>
struct scalar_power:
    public scalar_binary_functor<M, P> {
    typedef typename scalar_binary_functor<M, P>::argument1_type argument1_type;
    typedef typename scalar_binary_functor<M, P>::argument2_type argument2_type;
    typedef typename scalar_binary_functor<M, P>::result_type result_type;

    static BOOST_UBLAS_INLINE
    result_type apply (argument1_type t1, argument2_type t2) {
        return pow(abs(t1), t2);
    }
};

template<class M, class P>
BOOST_UBLAS_INLINE
typename enable_if< is_convertible<P, typename M::value_type>,
typename matrix_binary_scalar2_traits<M, const P, scalar_power<typename M::value_type, P> >::result_type
>::type
operator ^ (const matrix_expression<M> &m,
            const P &p) {
    typedef typename matrix_binary_scalar2_traits<M, const P, scalar_power<typename M::value_type, P> >::expression_type expression_type;
    return expression_type (m(), p);
}
}}}

After this your expression becomes:

D = element_prod(A^2, element_prod(B^3, C));